Xiaohua Luo
Junjie Guo
Zhongzhou Si
Jiequn Li- Department of Liver Transplant, The Second Xiangya Hospital, Central South University, Changsha, China
Backgrounds: Non-alcoholic fatty liver disease (NAFLD) presents as a common liver disease characterized by an indistinct pathogenesis. Disulfidptosis is a recently identified mode of cell death. This study aimed to investigate the potential role of disulfidptosis-related genes (DRGs) in the pathogenesis of NAFLD.
Methods: Gene expression profiles were obtained from the bulk RNA dataset GSE126848 and the single-cell RNA dataset GSE136103, both associated with NAFLD. Our study assessed the expression of DRGs in NAFLD and normal tissues. Weighted gene co-expression network analysis (WGCNA) and differential expression analysis were employed to identify the key NAFLD-specific differentially expressed DRGs (DE-DRGs). To explore the biological functions and immune regulatory roles of these key DE-DRGs, we conducted immune infiltration analysis, functional enrichment analysis, consensus clustering analysis, and single-cell differential state analysis. Finally, we validated the expression and biological functions of DRGs in NAFLD patients using histology and RNA-sequencing transcriptomic assays with human liver tissue samples.
Results: Through the intersection of WGCNA, differentially expressed genes, and DRGs, two key DE-DRGs (DSTN and MYL6) were identified. Immune infiltration analysis indicated a higher proportion of macrophages, T cells, and resting dendritic cells in NAFLD compared to control liver samples. Based on the key DE-DRGs, Two disulfidptosis clusters were defined in GSE126848. Cluster 1, with higher expression of the key DE-DRGs, exhibited increased immune infiltration abundance and was closely associated with oxidative stress and immune regulation compared to cluster 2. High-resolution analysis of mononuclear phagocytes highlighted the potential role of MYL6 in intrahepatic M1 phenotype Kupffer cells in NAFLD patients. Our transcriptome data revealed that the expression levels of the majority of DRGs were significantly increased in NAFLD patients. NAFLD patients exhibit elevated MYL6 correlating with inflammation, oxidative stress, and disease severity, offering promising diagnostic specificity.
Conclusion: This comprehensive study provides evidence for the association between NAFLD and disulfidptosis, identifying potential target genes and pathways in NAFLD. The identification of MYL6 as a possible treatment target for NAFLD provided a novel understanding of the disease’s development.
1 Introduction
Non-alcoholic fatty liver disease (NAFLD) comprises a collection of chronic liver disorders distinguished by the excessive accumulation of lipids within hepatocytes and the presence of steatosis (1). The global occurrence of this condition has significantly risen over the years, making NAFLD the prevailing liver ailment across the globe (2). It is currently estimated to affect up to 25% of the population (3). It covers a range of liver conditions that vary from simple steatosis (non-alcoholic fatty liver, NAFL) to non-alcoholic steatohepatitis (NASH). Increasing epidemiological evidence suggests that NAFLD is quickly becoming a primary cause of hepatocellular carcinoma (HCC) in numerous instances (4, 5). Despite the increasing prevalence of NAFLD over time, its underlying causes remain unclear and there is a lack of definitive treatment options (6). Therefore, it is critically important to identify potential pathogenic and therapeutic targets in NAFLD.
Recent research has identified a novel metabolism-related regulated cell death mechanism termed disulfidptosis (7). Xiaoguang Liu et al. discovered that when the intracellular reducing molecule nicotinamide adenine dinucleotide phosphate was depleted, leading to the accumulation of disulfide compounds such as cystine and triggering disulfide stress. This led to the creation of disulfide bonds between cytoskeletal proteins and the collapse of the actin filament network, ultimately causing disulfidptosis to occur (7). In this process, the cellular redox state played a crucial role. It was widely recognized that individuals with NAFLD demonstrate hepatic oxidative stress caused by compromised mitochondrial respiratory capacity and proton leakage (8). Oxidative stress in NAFLD can induce the activation of Kupffer cells, leading to their polarization towards the pro-inflammatory M1 phenotype. This polarization causes the discharge of many pro-inflammatory cytokines, promoting an inflammatory reaction that harms liver cells and speeds up the advancement of NAFLD (9, 10). Furthermore, previous studies have indicated a close association between sulfide metabolism and the development of NAFLD (11, 12). All of these findings indicated that disulfidptosis-related genes (DRGs) might carry a crucial component in the development and progression of NAFLD.
In this research, we employed weighted gene co-expression network analysis (WGCNA) and differential expression analysis to identify 312 differentially expressed genes (DEGs) strongly associated with NAFLD, elucidating their enriched biological pathways. We also investigated the immune microenvironment in NAFLD patients. Subsequently, we identified two key NAFLD-specific differentially expressed DRGs (DE-DRGs), MYL6 and DSTN, based on 23 DRGs (7). We characterized two subtypes of disulfidptosis associated with NAFLD and immune infiltration. Gene set enrichment analysis (GSEA) results showed that oxidative stress pathways, triglyceride metabolism pathways, and pro-inflammatory pathways were activated in the subtype with high expression of key DE-DRGs. Furthermore, MYL6 was found to be highly expressed in liver immune cells and might potentially promote Kupffer cell polarization towards the pro-inflammatory M1 phenotype through oxidative stress. To reinforce our findings, we further validated the expression of the 23 DRGs in NAFLD patients and explored the associated biological functions and clinical significance of MYL6.
2 Materials and methods
2.1 Gene expression profiles collection
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.
Figure 1 Research design flow chart.
2.2 Weighted gene co-expression network analysis
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.
2.3 Differential expression analysis
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.
2.4 Functional enrichment analysis
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.
2.5 Immune infiltration analysis
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).
2.6 Consensus clustering analysis for the key DE-DRGs
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.
2.7 Analysis of scRNA-seq data
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).
2.8 Clinical samples and patient data
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.
2.9 Hematoxylin and eosin staining of liver tissue
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.
2.10 Liver transcriptome sequencing
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.
2.11 Statistical analysis
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.
3 Results
3.1 Building co-expression networks and screening of critical module
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 R2 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).
Figure 2 WGCNA analysis and clinically significant module identification. (A) Sample clustering dendrogram for outlier detection. (B) Clustering dendrogram of 57 liver samples with trait heatmap. (C) Under different soft thresholds, the scale-free fit index and mean connectivity are shown. (D) Hierarchical clustering of genes in the dendrogram. (E) The correlation between clinical traits and gene modules is represented as a heatmap. (F) A scatter plot of the yellow module’s eigengenes. (G) The co-expression modules are plotted as a network heatmap. (H) The eigengene network shows the relationships between modules and clinical trait weights.
3.2 Identification of NAFLD-specific DEGs and their KEGG enrichment analysis
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.
Figure 3 Identification of NAFLD-specific DEGs and their KEGG enrichment analysis. (A) Volcano plot of DEGs for the GSE126848 dataset. The top 5 most significantly up- and down-regulated DEGs are labeled individually. (B) Venn diagram showing the NAFLD-specific DEGs common to DEGs and WGCNA. (C) Heatmap of NAFLD-specific DEGs expression data. (D) Bar chart of KEGG enrichment analysis of NAFLD-specific DEGs. (E) Network diagram of pathway enrichment analysis.
3.3 The immune microenvironment in patients with 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.
Figure 4 Identification of the immune microenvironment in patients with NAFLD by CIBERSORT algorithm. (A) The relative ratio of immunocytes in patients with NAFLD and controls. (B) Proportional correlation heatmap of immune cells. The gradient color represents the correlation strength. (C) A box plot comparing the infiltration rates of 22 immune cells between NAFLD patients and controls. (D) Heatmap of the relative proportion of immune cells in the samples. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.
3.4 Identification of key DE-DRGs in NAFLD and disulfidptosis clusters immune infiltration analysis
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.
Figure 5 Identification of the key DE-DRGs in NAFLD. (A) Chromosome position of the 23 DRGs. (B) Venn diagram showing the overlapping genes common to DRGs, DEGs, and WGCNA. (C) Consensus matrix plots of the key DE-DRGs. (D) Delta area plot. (E) Heatmap of relative expression distribution of key DE-DRGs in two clusters. (F) Box plot showing immune infiltration differences between disulfidptosis clusters by CIBERSORT algorithm. (G) Correlation analysis of the key DE-DRGs with immune cells. (H) GSEA analysis based on the canonical pathways gene sets. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.
3.5 Processing of the GSE136103 single-cell data set and construction of a cell atlas for NAFLD
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).
Figure 6 Clustering of GSE136103 scRNA-seq data and identification of cell types. (A) Violin plot of quality control parameters for scRNA-seq of liver tissue samples. (B) Heatmap of each cell cluster’s top 60 specific marker genes, while the upper panel of the heatmap shows the hierarchical clustering of clusters and the number of cells. In the heatmap, the gene expression is presented after Z-score processing, and the redder color represents the higher gene expression. (C) A bar chart shows the cell composition of cell clusters in 8 liver tissue samples. (D) UMAP plot (center) shows 23,929 cells from 8 samples, colored by major classes. The circle is divided into three parts from inside to outside, representing the distribution of samples in each cell cluster, the distribution of each cell cluster between groups, and the name of the cell clusters.
3.6 The key DE-DRGs displayed significant differential expression in the MPs cluster
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.
Figure 7 Weighted kernel density estimation of the key DE-DRGs and differential state analysis in GSE136103. (A) UMAP plot showing the Nebulosa expression densities for MYL6. (B) UMAP plot showing the Nebulosa expression densities for DSTN. (C) Pseudobulk-level multi-dimensional scaling (MDS) plot. (D) The violin plot shows the expression of MYL6 for the mononuclear phagocyte cluster. (E) The volcano plot shows the results of the differential state analysis in GSE136103 with the key DE-DRGs labeled.
3.7 MYL6 might be involved in Kupffer cell polarization in NAFLD patients
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.
Figure 8 The specific phenotypes of mononuclear phagocytes. (A) UMAP plot showing 5 distinct clusters in mononuclear phagocytes, colored by cluster ID (left panel). A distance matrix generated in either gene expression space is used to infer the phylogenetic tree of 5 clusters. (right panel). (B) UMAP plot showing the Nebulosa expression densities for MYL6 in mononuclear phagocytes. (C) Heatmap of marker genes among different clusters. (D) Violin plots showing the functional module (Mod) scores of mononuclear phagocyte subsets. (E) Pseudobulk-level MDS plot. (F) The volcano plot shows the results of the differential state analysis in each mononuclear phagocyte cluster with MYL6 labeled. (G) GSEA in C02-CD5L-Macro cluster based on the canonical pathways gene sets.
3.8 The expression of MYL6 in patients with NAFLD and its close correlation with clinical parameters
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).
Figure 9 Validation of MYL6 in our data. (A) Liver tissue sections stained with H&E. (B) Relative expression of DRGs in liver samples from patients with NAFLD and controls. (C) Heatmap showing the MYL6 was significantly positively correlated with the expression of Kupffer cell signature genes, M1 phenotype macrophage signature genes, and oxidative stress-related genes. (D) GSEA in NAFLD patients based on the median MYL6 expression. (E, F) Box plot showing the expression of MYL6 among the control, NAFL, and NASH groups based on the different data. (G) Relative expression of MYL6 in liver tissue samples from patients with NAFLD-LIHC and controls. (H) Validation of MYL6 in NAFLD diagnosis in our data with ROC curve analysis. (I–M) Correlations of intrahepatic MYL6 mRNA levels with TG, TC, LDL-C, HDL-C, and NAS in our data. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. ns: no statistically significant (P > 0.05).
4 Discussion
Due to the escalating prevalence of obesity and weight-related metabolic complications, NAFLD has become the most prevalent chronic liver ailment worldwide (41). In a subset of patients, hepatic steatosis progresses to NASH and may further develop into fibrosis, cirrhosis, and NAFLD-related HCC (42–44). Fatty acid metabolic reprogramming unique to hepatocytes is a key indicator of the liver cancer (45). Despite extensive research, the pathogenesis of NAFLD remains complex and not entirely understood, leading to a lack of effective therapeutic drugs, reliable non-invasive diagnostic tools, and dynamic biomarkers (46). Consequently, the identification of pathogenic targets and effective biomarkers for NAFLD is crucial for formulating personalized treatment strategies. There is no consensus on the role of disulfideptosis in the occurrence of NAFLD and its progression as a novel form of cell death (7). We aimed to investigate the potential role of DRGs in the pathogenesis of NAFLD and immune microenvironment based on single-cell and bulk RNA-seq.
In our study, a yellow module with 564 genes was identified, which has the highest correlation with the phenotype (NAFLD) in WGCNA. Furthermore, we identified 312 NAFLD-specific DEGs in DEGs and WGCNA, of which 293 were up-regulated and 19 were down-regulated. Functional enrichment analysis in the present study found that the NAFLD-specific DEGs were associated with important pathogenesis processes and pathways of NAFLD (39, 40), including Non−alcoholic fatty liver disease pathway, Oxidative phosphorylation pathway, and PPAR signaling pathways. These findings strongly indicated that the NAFLD-specific DEGs identified in this study play an important role in the occurrence and development of NAFLD and should be further investigated. Subsequently, immune infiltration analysis revealed significant alterations in the composition of immune cells in patients with NAFLD. The percentage of plasma cells, CD8 T cells, gamma delta T cells, M1 macrophages, and resting dendritic cells was higher in NAFLD samples compared to control samples, while a lower percentage of M2 macrophages was present. These were consistent with the results of previous studies (47–51). Although NAFLD is primarily a metabolic disorder, it also involves immune cell-mediated inflammatory processes. Inflammation becomes an indispensable factor in disease progression, particularly in the NASH stage. Liver immune cells are diverse in their steady state and undergo further evolution during NAFLD, directly impacting the severity of the disease (48).
To explore the role of disulfidptosis features in NAFLD, we examined the expression patterns of DRGs in the control and NAFLD groups. The results revealed that most of the DRGs were significantly differentially expressed in patients with NAFLD, which suggested a crucial role for DRGs in the occurrence and progression of NAFLD. On the other hand, we employed unsupervised clustering to categorize patients into two distinct disulfidptosis subtypes. We found that the key NAFLD-specific DE-DRGs (DSTN and MYL6) were up-regulated in cluster 1. Interestingly, cluster 1 has a high abundance of immune cell infiltration. GSEA analysis indicated that in cluster 1, oxidative stress-related pathways and pro-inflammatory pathways were upregulated. The expression levels of key DE-DRGs were positively correlated with the abundance of M1 macrophages and CD8 T cell infiltration. Prior studies have also found that disulfidptosis-related long non-coding RNAs are closely associated with the level of immune cell infiltration in liver tissue (52). Additionally, previous studies showed that elevated levels of reactive oxygen species played a crucial role in the inflammatory response, fibrosis, necrosis, and apoptosis occurring in NAFLD (53). Oxidative stress can induce Kupffer cells to produce various cytokines, such as TNF-α, thereby exacerbating inflammation and cell apoptosis (10). As for HSCs, their proliferation and collagen synthesis are triggered by lipid peroxidation induced by oxidative stress (54). This suggested that DRGs may promote immune cell infiltration in NAFLD through oxidative stress-related pathways, thereby accelerating the progression of the disease.
Subsequently, we further explored the expression distribution and biological functions of key DE-DRGs at the single-cell level within the liver. Traditional transcriptomic data from bulk RNA-seq pose challenges in depicting the heterogeneity of various cell types within the livers of both NAFLD patients and healthy individuals. With technological advancements, high-throughput sequencing techniques, such as scRNA-seq, have surfaced to furnish transcriptomic information at the cellular level (55, 56). We annotated and identified ten different cell subtypes based on scRNA-seq data. Results revealed that MYL6 was expressed in the majority of immune cells within the liver, while DSTN was predominantly expressed in cholangiocytes and ECs. Furthermore, MYL6 exhibited significant differential expression in the MPs cluster, with its expression trend aligning with that observed in bulk RNA samples. Subsequently, further refinement of the MPs cell subtype revealed that MYL6 was predominantly expressed in the CD5L-Macro resident in the livers of NAFLD patients. Functional module scores and enrichment analysis unveiled that the CD5L-Macro was characterized by high phagocytic activity and secretion of pro-inflammatory and fibrotic mediators, belonging to the M1 phenotype kupfer cells. GSEA enrichment analysis of C02-CD5L-Macro indicated activation of pro-inflammatory and pro-fibrotic signaling pathways, pathways related to regulatory cell death (such as regulated necrosis and ferroptosis), lipid metabolism pathways, and pathways associated with sulfur metabolism (such as sulfate formation due to sulfide oxidation) in NAFLD patients compared to the controls. Numerous studies suggested that Kupffer cells, a self-sustaining macrophage population residing in the liver, played a role in the progression of NAFLD (57, 58). During the onset of NAFLD, Kupffer cells underwent pro-inflammatory polarization, serving as crucial inflammatory mediators (51, 59). Subsequently, cytokines released by Kupffer cells stimulate a pro-inflammatory cascade, activating other immune cells (such as T cells), inducing the expression of transcription factors involved in lipid metabolism and transport, and promoting programmed cell death (60, 61). Therefore, we have reason to believe that MYL6 may contribute to the development of NAFLD through the promotion of M1 phenotype Kupffer cells.
Finally, we designed a validation study to investigate whether DRGs exhibited differential expression in NAFLD across 16 human liver tissue samples. Based on the RNA-seq results, the expression of the majority of DRGs was significantly increased in the NAFLD patients compared to the controls. In this study, particular emphasis was placed on the investigation of MYL6. Limited knowledge exists regarding the role of MYL6 in NAFLD. MYL6 encodes a myosin alkali light chain that is expressed in smooth muscle and non-muscle tissues, which is important for the assembly of cytoskeletal structures (62). Previous studies have identified a close association between MYL6 and the inflammatory disease Multiple Sclerosis, suggesting its predictive value for prognosis (63). Furthermore, MYL6 exhibited elevated expression in conditions such as obesity, asthma, and cervical cancer, yet the underlying mechanisms remain not fully elucidated (64). Our transcriptomic analysis revealed upregulation of MYL6 expression in the NAFLD group, closely associated with the expression of characteristic genes of Kupffer cells, M1 phenotype macrophages, and oxidative stress-related genes. The high-expression group of MYL6 exhibited upregulation of pathways related to inflammation, fatty acid metabolism, and oxidative stress. Intriguingly, the expression of MYL6 was significantly positively correlated with the progression of the disease spectrum in NAFLD. Evaluation of clinical features indicated a significant positive correlation between MYL6 expression levels and NAS, as well as concentrations of TC, TG, and LDL-C in the serum, while a significant negative correlation was observed with serum HDL-C concentration. Our findings suggested that MYL6 may promote the occurrence of oxidative stress, increase the infiltration of M1 phenotype Kupffer cells, and facilitate the progression of NAFLD through multiple pathways.
Certainly, this study inevitably has some limitations. The investigation solely validated the mRNA levels of DRGs and did not explore their protein-level changes in NAFLD. Additionally, further experimental validation was required to elucidate the potential regulatory role of DRGs in the occurrence and development of NAFLD. Finally, larger sample sizes in subsequent studies were needed to clarify the correlation between MYL6 expression levels and clinical pathological indicators of NAFLD.
In summary, this study comprehensively demonstrated the relationship between NAFLD and disulfidptosis using both single-cell and bulk RNA-seq data. It elucidated the regulatory role of DRGs in the hepatic immune microenvironment of NAFLD patients. Numerous potential target genes and pathways have been identified in NAFLD, providing potential targets for therapeutic interventions. Importantly, we revealed MYL6 as a novel immunomodulator participating in the pathogenesis of NAFLD. In conclusion, our research findings may offer a better understanding of the underlying mechanisms, thereby providing novel targets for the treatment of NAFLD.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GSE260666 (GEO).
Ethics statement
The studies involving humans were approved by the Ethics Committee of the Second Xiangya Hospital of Central South University (No. 2019-050). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.
Author contributions
XL: Writing – original draft, Validation, Methodology, Investigation, Formal analysis, Data curation, Visualization. JG: Writing – original draft, Methodology, Investigation. HD: Writing – original draft, Resources, Data curation. ZH: Visualization, Software, Writing – original draft, Methodology. YW: Methodology, Investigation, Writing – original draft, Resources. ZS: Writing – review & editing, Supervision, Funding acquisition. JL: Writing – review & editing, Visualization, Supervision, Project administration, Funding acquisition, Conceptualization.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded by the Key Research and Development Program of Hunan (Grant No. 2022SK2040), National Natural Science Foundation of China (Grant No. 82070679), and Natural Science Foundation of Hunan province (Grant No. 2022JJ30813 and 2020JJ4804).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1386905/full#supplementary-material
References
1. Brunt EM. Pathology of nonalcoholic fatty liver disease. Nat Rev Gastroenterol Hepatol. (2010) 7:195–203. doi: 10.1038/nrgastro.2010.21
2. Lazarus JV, Ekstedt M, Marchesini G, Mullen J, Novak K, Pericas JM, et al. A cross-sectional study of the public health response to non-alcoholic fatty liver disease in Europe. J Hepatol. (2020) 72:14–24. doi: 10.1016/j.jhep.2019.08.027
3. Younossi ZM, Koenig AB, Abdelatif D, Fazel Y, Henry L, Wymer M. Global epidemiology of nonalcoholic fatty liver disease-meta-analytic assessment of prevalence, incidence, and outcomes. Hepatology. (2016) 64:73–84. doi: 10.1002/hep.28431
4. Anstee QM, Reeves HL, Kotsiliti E, Govaere O, Heikenwalder M. From Nash to Hcc: current concepts and future challenges. Nat Rev Gastroenterol Hepatol. (2019) 16:411–28. doi: 10.1038/s41575–019-0145–7
5. Sanyal AJ. Past, present and future perspectives in nonalcoholic fatty liver disease. Nat Rev Gastroenterol Hepatol. (2019) 16:377–86. doi: 10.1038/s41575–019-0144–8
6. Seebacher F, Zeigerer A, Kory N, Krahmer N. Hepatic lipid droplet homeostasis and fatty liver disease. Semin Cell Dev Biol. (2020) 108:72–81. doi: 10.1016/j.semcdb.2020.04.011
7. Liu X, Nie L, Zhang Y, Yan Y, Wang C, Colic M, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol. (2023) 25:404–14. doi: 10.1038/s41556–023-01091–2
8. Koliaki C, Szendroedi J, Kaul K, Jelenik T, Nowotny P, Jankowiak F, et al. Adaptation of hepatic mitochondrial function in humans with non-alcoholic fatty liver is lost in steatohepatitis. Cell Metab. (2015) 21:739–46. doi: 10.1016/j.cmet.2015.04.004
9. Tanaka M, Miyajima A. Liver regeneration and fibrosis after inflammation. Inflammation Regener. (2016) 36:19. doi: 10.1186/s41232–016-0025–2
10. Tang S-P, Mao X-L, Chen Y-H, Yan L-L, Ye L-P, Li S-W. Reactive oxygen species induce fatty liver and ischemia-reperfusion injury by promoting inflammation and cell death. Front In Immunol. (2022) 13:870239. doi: 10.3389/fimmu.2022.870239
11. Chen M, Zhu J-Y, Mu W-J, Luo H-Y, Li Y, Li S, et al. Cdo1-Camkk2-Ampk axis confers the protective effects of exercise against Nafld in mice. Nat Commun. (2023) 14:8391. doi: 10.1038/s41467–023-44242–7
12. Dhamija RK, Gaba P, Arora S, Kaintura A, Kumar M, Bhattacharjee J. Homocysteine and lipoprotein (a) correlation in ischemic stroke patients. J Neurol Sci. (2009) 281:64–8. doi: 10.1016/j.jns.2009.02.341
13. Suppli MP, Rigbolt KTG, Veidal SS, Heebøll S, Eriksen PL, Demant M, et al. Hepatic transcriptome signatures in patients with varying degrees of nonalcoholic fatty liver disease compared with healthy normal-weight individuals. Am J Physiol Gastrointest Liver Physiol. (2019) 316:G462–G72. doi: 10.1152/ajpgi.00358.2018
14. Ramachandran P, Dobie R, Wilson-Kanamori JR, Dora EF, Henderson BEP, Luu NT, et al. Resolving the fibrotic niche of human liver cirrhosis at single-cell level. Nature. (2019) 575:512–8. doi: 10.1038/s41586–019-1631–3
15. Langfelder P, Horvath S. Wgcna: an R package for weighted correlation network analysis. BMC Bioinf. (2008) 9:559. doi: 10.1186/1471–2105-9–559
16. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. (2010) 11:R106. doi: 10.1186/gb-2010–11-10-r106
17. Kanehisa M, Goto S, Kawashima S, Okuno Y, Hattori M. The Kegg resource for deciphering the genome. Nucleic Acids Res. (2004) 32:D277–80. doi: 10.1093/nar/gkh063
18. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). (2021) 2:100141. doi: 10.1016/j.xinn.2021.100141
19. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database (Msigdb) hallmark gene set collection. Cell Syst. (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004
20. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118
21. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. (2019) 37:773–82. doi: 10.1038/s41587–019-0114–2
22. Hänzelmann S, Castelo R, Guinney J. Gsva: gene set variation analysis for microarray and Rna-Seq data. BMC Bioinf. (2013) 14:7. doi: 10.1186/1471–2105-14–7
23. Aran D, Hu Z, Butte AJ. Xcell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059–017-1349–1
24. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, et al. Comprehensive integration of single-cell data. Cell. (2019) 177:1888–1902.e21. doi: 10.1016/j.cell.2019.05.031
25. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with harmony. Nat Methods. (2019) 16:1289–96. doi: 10.1038/s41592–019-0619–0
26. Pallett LJ, Davies J, Colbeck EJ, Robertson F, Hansi N, Easom NJW, et al. Il-2high tissue-resident T cells in the human liver: sentinels for hepatotropic infection. J Exp Med. (2017) 214:1567–80. doi: 10.1084/jem.20162115
27. Guo X, Zhang Y, Zheng L, Zheng C, Song J, Zhang Q, et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nat Med. (2018) 24:978–85. doi: 10.1038/s41591–018-0045–3
28. Zhang C, Li J, Cheng Y, Meng F, Song J-W, Fan X, et al. Single-cell Rna sequencing reveals intrahepatic and peripheral immune characteristics related to disease phases in Hbv-infected patients. Gut. (2023) 72:153–67. doi: 10.1136/gutjnl-2021–325915
29. Zhang X, Lan Y, Xu J, Quan F, Zhao E, Deng C, et al. Cellmarker: A manually curated resource of cell markers in human and mouse. Nucleic Acids Res. (2019) 47:D721–D8. doi: 10.1093/nar/gky900
30. Crowell HL, Soneson C, Germain P-L, Calini D, Collin L, Raposo C, et al. Muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data. Nat Commun. (2020) 11:6077. doi: 10.1038/s41467–020-19894–4
31. Bischoff P, Trinks A, Obermayer B, Pett JP, Wiederspahn J, Uhlitz F, et al. Single-cell Rna sequencing reveals distinct tumor microenvironmental patterns in lung adenocarcinoma. Oncogene. (2021) 40:6748–58. doi: 10.1038/s41388–021-02054–3
32. Li X, Li S, Wu B, Xu Q, Teng D, Yang T, et al. Landscape of immune cells heterogeneity in liver transplantation by single-cell Rna sequencing analysis. Front In Immunol. (2022) 13:890019. doi: 10.3389/fimmu.2022.890019
33. Guillot A, Tacke F. Liver macrophages: old dogmas and new insights. Hepatol Commun. (2019) 3:730–43. doi: 10.1002/hep4.1356
34. Aizarani N, Saviano A, Sagar, Mailly L, Durand S, Herman JS, et al. A human liver cell atlas reveals heterogeneity and epithelial progenitors. Nature. (2019) 572:199–204. doi: 10.1038/s41586–019-1373–2
35. Wu H, Gonzalez Villalobos R, Yao X, Reilly D, Chen T, Rankin M, et al. Mapping the single-cell transcriptomic response of murine diabetic kidney disease to therapies. Cell Metab. (2022) 34:1064–78.e6. doi: 10.1016/j.cmet.2022.05.010
36. Alquicira-Hernandez J, Powell JE. Nebulosa recovers single-cell gene expression signals by kernel density estimation. Bioinformatics. (2021) 37:2485–7. doi: 10.1093/bioinformatics/btab003
37. La Manno G, Siletti K, Furlan A, Gyllborg D, Vinsland E, Mossi Albiach A, et al. Molecular architecture of the developing mouse brain. Nature. (2021) 596:92–6. doi: 10.1038/s41586-021-03775-x
38. Kleiner DE, Brunt EM, Van Natta M, Behling C, Contos MJ, Cummings OW, et al. Design and validation of a histological scoring system for nonalcoholic fatty liver disease. Hepatology. (2005) 41:1313–21. doi: 10.1002/hep.20701
39. Loomba R, Friedman SL, Shulman GI. Mechanisms and disease consequences of nonalcoholic fatty liver disease. Cell. (2021) 184:2537–64. doi: 10.1016/j.cell.2021.04.015
40. Bessone F, Razori MV, Roma MG. Molecular pathways of nonalcoholic fatty liver disease development and progression. Cell Mol Life Sci. (2019) 76:99–128. doi: 10.1007/s00018–018-2947–0
41. Angulo P. Nonalcoholic fatty liver disease. New Engl J Med. (2002) 346:1221–31. doi: 10.1056/NEJMra011775
42. Llovet JM, Willoughby CE, Singal AG, Greten TF, Heikenwälder M, El-Serag HB, et al. Nonalcoholic steatohepatitis-related hepatocellular carcinoma: pathogenesis and treatment. Nat Rev Gastroenterol Hepatol. (2023) 20:487–503. doi: 10.1038/s41575–023-00754–7
43. Xu K, Wu T, Xia P, Chen X, Yuan Y. Alternative splicing: A bridge connecting Nafld and Hcc. Trends Mol Med. (2023) 29:859–72. doi: 10.1016/j.molmed.2023.07.001
44. Xu K, Xia P, Liu P, Zhang X. A six lipid metabolism related gene signature for predicting the prognosis of hepatocellular carcinoma. Sci Rep. (2022) 12:20781. doi: 10.1038/s41598–022-25356–2
45. Xu K, Xia P, Chen X, Ma W, Yuan Y. Ncrna-mediated fatty acid metabolism reprogramming in Hcc. Trends In Endocrinol Metab: TEM. (2023) 34:278–91. doi: 10.1016/j.tem.2023.02.007
46. Wu T, Chen X, Xu K, Dai C, Li X, Zhang Y-W-Q, et al. Lim domain only 7 negatively controls nonalcoholic steatohepatitis in the setting of hyperlipidemia. Hepatol (Baltimore Md). (2024) 79:149–66. doi: 10.1097/HEP.0000000000000585
47. Deczkowska A, David E, Ramadori P, Pfister D, Safran M, Li B, et al. Xcr1+ Type 1 conventional dendritic cells drive liver pathology in non-alcoholic steatohepatitis. Nat Med. (2021) 27:1043–54. doi: 10.1038/s41591–021-01344–3
48. Huby T, Gautier EL. Immune cell-mediated features of non-alcoholic steatohepatitis. Nat Rev Immunol. (2022) 22:429–43. doi: 10.1038/s41577–021-00639–3
49. Yang T, Qu X, Wang X, Xu D, Sheng M, Lin Y, et al. Macrophage sting-yap axis controls hepatic steatosis by promoting autophagic degradation of lipid droplets. Hepatology. (2023). doi: 10.1097/HEP.0000000000000638
50. Shalapour S, Lin X-J, Bastian IN, Brain J, Burt AD, Aksenov AA, et al. Inflammation-induced Iga+ Cells dismantle anti-liver cancer immunity. Nature. (2017) 551:340–5. doi: 10.1038/nature24302
51. Krenkel O, Tacke F. Liver macrophages in tissue homeostasis and disease. Nat Rev Immunol. (2017) 17:306–21. doi: 10.1038/nri.2017.11
52. Xu K, Dai C, Yang J, Xu J, Xia C, Li J, et al. Disulfidptosis-related Lncrna signatures assess immune microenvironment and drug sensitivity in hepatocellular carcinoma. Comput Biol Med. (2024) 169:107930. doi: 10.1016/j.compbiomed.2024.107930
53. Li S, Tan H-Y, Wang N, Zhang Z-J, Lao L, Wong C-W, et al. The role of oxidative stress and antioxidants in liver diseases. Int J Mol Sci. (2015) 16:26087–124. doi: 10.3390/ijms161125942
54. Allameh A, Niayesh-Mehr R, Aliarab A, Sebastiani G, Pantopoulos K. Oxidative stress in liver pathophysiology and disease. Antioxid (Basel). (2023) 12:1653. doi: 10.3390/antiox12091653
55. Papalexi E, Satija R. Single-cell Rna sequencing to explore immune cell heterogeneity. Nat Rev Immunol. (2018) 18:35–45. doi: 10.1038/nri.2017.76
56. Van de Sande B, Lee JS, Mutasa-Gottgens E, Naughton B, Bacon W, Manning J, et al. Applications of single-cell Rna sequencing in drug discovery and development. Nat Rev Drug Discovery. (2023) 22:496–520. doi: 10.1038/s41573–023-00688–4
57. Zhang X, Han J, Man K, Li X, Du J, Chu ESH, et al. Cxc chemokine receptor 3 promotes steatohepatitis in mice through mediating inflammatory cytokines, macrophages and autophagy. J Hepatol. (2016) 64:160–70. doi: 10.1016/j.jhep.2015.09.005
58. Kazankov K, Barrera F, Møller HJ, Rosso C, Bugianesi E, David E, et al. The macrophage activation marker Scd163 is associated with morphological disease stages in patients with non-alcoholic fatty liver disease. Liver Int. (2016) 36:1549–57. doi: 10.1111/liv.13150
59. Huang W, Metlakunta A, Dedousis N, Zhang P, Sipula I, Dube JJ, et al. Depletion of liver Kupffer cells prevents the development of diet-induced hepatic steatosis and insulin resistance. Diabetes. (2010) 59:347–57. doi: 10.2337/db09–0016
60. Vonderlin J, Chavakis T, Sieweke M, Tacke F. The multifaceted roles of macrophages in Nafld pathogenesis. Cell Mol Gastroenterol Hepatol. (2023) 15:1311–24. doi: 10.1016/j.jcmgh.2023.03.002
61. Zhang W, Lang R. Macrophage metabolism in nonalcoholic fatty liver disease. Front In Immunol. (2023) 14:1257596. doi: 10.3389/fimmu.2023.1257596
62. Vierthaler M, Sun Q, Wang Y, Steinfass T, Poelchen J, Hielscher T, et al. Adck2 knockdown affects the migration of melanoma cells via Myl6. Cancers (Basel). (2022) 14:1071. doi: 10.3390/cancers14041071
63. Ye F, Liang J, Li J, Li H, Sheng W. Development and validation of a five-gene signature to predict relapse-free survival in multiple sclerosis. Front Neurol. (2020) 11:579683. doi: 10.3389/fneur.2020.579683
Keywords: NAFLD, disulfidptosis, oxidative stress, Kupffer cell, MYL6
Citation: Luo X, Guo J, Deng H, He Z, Wen Y, Si Z and Li J (2024) Unveiling the role of disulfidptosis-related genes in the pathogenesis of non-alcoholic fatty liver disease. Front. Immunol. 15:1386905. doi: 10.3389/fimmu.2024.1386905
Received: 16 February 2024; Accepted: 03 May 2024;
Published: 15 May 2024.
Edited by:
Jia Xiao, Jinan University, ChinaCopyright © 2024 Luo, Guo, Deng, He, Wen, Si and Li. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhongzhou Si, zhongzsi@csu.edu.cn; Jiequn Li, leejiequn@csu.edu.cn