Exploring the biological functions and immune regulatory roles of IRAK3, TNFRSF1A, CX3CR1, and JUNB in T2DM combined with MAFLD: integrated bioinformatics and single-cell analysis

The datasets GSE15653, GSE89632, GSE24807 and GSE23343 were obtained from the GEO database through the GEOquery package (version 2.72.0) (6),please refer to Supplementary Table S1 for more details.GSE15653 and GSE89632 served as training sets, while GSE24807 and GSE23343 was utilized as the validation set. The MAFLD samples within the GSE89632 dataset were categorized by different stages: Simple Steatosis (SS) and Non-Alcoholic Steatohepatitis (NASH), and were designated as GSE89632_Simple Steatosis (GSE89632_SS) and GSE89632_Non-Alcoholic Steatohepatitis (GSE89632_NASH), respectively. The data from GSE15653, GSE89632_SS, and GSE89632_NASH were normalized using the normalizeBetweenArrays function from the limma package (7). Principal component analysis (PCA) was performed on the normalized dataset. Difference analysis was performed separately to obtain the respective DEGs, using limma package (version 3.58.1), with screening criteria: |log2FC| > 0.5 and P < 0.05 (8, 9). Subsequent analyses included the generation of PCA plots, volcano plots, and basic numerical heatmaps using "ggplot2," alongside the creation of complex numerical heatmaps utilizing the "ComplexHeatmap" package (version 2.20.0) (10).

Weighted gene co-expression networks were analyzed using the WGCNA package (version 1.72-5) (11) for the datasets GSE89632_SS and GSE89632_NASH. Initially, the analysis focused on the top 25% of genes based on expression variance, identifying and excluding potential outlier samples through clustering analysis. Next, the pickSoftThreshold function was employed to ascertain the optimal soft threshold. After establishing the threshold, the cutreeDynamic function identified dynamic modules and set the minimum gene count required for each module (12). Initial modules were identified through the Dynamic Tree Cut (DTC) method, with similar modules subsequently merged based on the clustering relationships of modular eigengene (ME). Subsequently, the TOMsimilarity function calculated the topological molecular similarity matrix among genes, with 1,000 genes randomly selected for heatmap visualization. Pearson correlation coefficients were computed to evaluate associations between modules and clinical traits, with statistical significance assessed via the corPvalueStudent function (WGCNA R package). Modules demonstrating significant correlations with both SS and NASH (p < 0.05) were selected for downstream analysis. KEGG pathway enrichment analysis was then systematically performed on the genes within these modules.

The results of DEGs of GSE15653 were intersected with DEGs of GSE89632_SS, GSE89632_NASH, and WGCNA modular genes, respectively, and the results were visualized in an "UpSet plot". In order to explore the potential interactions of the above genes, we used the STRING database (https://string-db.org/↗) for the analysis (13), and the results were imported into Cytoscape (version 3.10.2) for the analysis of protein interaction networks (14). Meanwhile, the Maximal Clique Centrality (MCC), Degree and Molecular Network Centrality (MNC) algorithms in CytoHubba plug-in were applied to select the top 15 genes at key positions in the PPI network (15). Subsequently, the results of the three algorithms were intersected and visualized using a Wayne diagram. Chromosomal localization analysis of the intersection hub genes was performed using the "circlize package". Finally, the correlation analysis was performed using Spearman, and the correlation pie charts were used to show the degree of association and interactions between genes, which were visualized using the "ggplot2 package".

In order to deeply investigate the molecular mechanisms related to autophagy in T2DM combined with MAFLD, we combined the resources from multiple databases. We obtained the autophagy-associated molecular mechanisms from PubMedGene database (https://www.ncbi.nlm.nih.gov/gene/↗), Genecard database (https://www.genecards.org/↗), and GSEA database (https://www.gsea-msigdb.org/↗), autophagy-related websites (http://www.autophagy.lu/↗), (http://hamdb.scbdd.com/home/index/↗) to obtain autophagy-related genes.

The GSE89632_SS and GSE89632NASH datasets underwent individual screenings using LASSO, RFE-SVM, and Boruta methods, respectively. The diagnostic LASSO coefficient screening leveraged the glmnet package (Version 4.1.7) to analyze the cleaned data, derive variable lambda, likelihood, and L1 regularization values, along with visualization. We employed ten-fold cross-validation (seed number: 2022) for validation (16). Additionally, we utilized two feature selection strategies: Recursive Feature Elimination (RFE) and the Boruta algorithm. RFE eliminates features without significant impact on Accuracy during model training repetitions (17), while the Boruta algorithm selects or eliminates significant and non-significant features in each base learner training iteration based on shadow features and Z-scores (18). This study integrated SVM base learner scores from the e1071 package (Version 1.7-14) with the RFE strategy (17) through 10*10 fold cross-validation and applied the Boruta algorithm using the Boruta package (Version 8.0.0) (18). Subsequently, the outcomes of Lasso, Boruta, and RFE-SVM screenings, combined with autophagy-related genes, were intersected to identify diagnostic biomarkers linked to autophagy in T2DM and NAFLD. These results were visually represented in Wayne diagrams.

In datasets GSE89632_SS, GSE89632_NASH and GSE15653, ROC curve analysis was performed on each of the three datasets using the pROC software package to determine the sensitivity and specificity of the above genes (19), to predict the ROC-related information and data of the variables at their respective cut-off values, and to assess the diagnostic accuracy, and the results were presented as ROC The results were quantified by area under the curve (AUC), and the AUC ranges between 0.5 and 1, and values closer to 1 indicate superior diagnostic performance of the variable in predicting clinical outcomes and visualized using ggplot2.

To investigate the functions and pathways of the above biomarkers, we used the DAVID online database (https://david.ncifcrf.gov/↗) for GO and KEGG enrichment analysis (20 –22). The "org.Hs.eg.db" package was used to convert the ID of the input gene list, and the "clusterProfiler package" was used for enrichment analysis (23). The significance of the enrichment results was assessed by calculating the z-score value for each enriched entry using the "GOplot" package (24). The results were screened for significance and biological significance at P<0.05 and FDR<0.2 (25). At the same time, GSEA(Gene Set Enrichment Analysis) was performed (26), and the results were screened according to the following criteria: normalized enrichment score |NES| > 1, FDR < 0.25, p.adjust < 0.05. Finally, the results were visualized using the "ggplot2 package".

In order to explore the immune cell infiltration in liver tissues of patients with T2DM combined with MAFLD, the "Single Sample Gene Set Enrichment Analysis" (ssGSEA) method was used in this study (27). The ssGSEA algorithm provided by the GSVA package was utilized for this purpose. Using markers specific to each class of immune cells (28) as gene sets, the enrichment score for each class of immune cells in each sample was calculated to assess the infiltration of immune cells in each sample. All analyses and visualizations were performed in R 4.2.1. We used the ggplot2 package to draw histograms to visualize the differences in immune cell infiltration between the normal and control groups. In addition, Spearman statistics were used to analyze the correlation between them one by one, and the linkET package was used for the calculation of the data portion of the network graph. The analysis results were visualized by ggplot2 package for group comparison plots, lollipop plots, correlation scatter plots, and correlation network heatmaps, thus demonstrating more intuitively the infiltration of immune cells in liver tissues of patients with T2DM combined with MAFLD.

Thirty healthy clean-grade male Sprague–Dawley rats, aged 6–8 weeks, weighing 180–220 g, were provided by Wuhan Yunkron Technology Company. The animal production and use license number is DCXR(E)2018-0021. SYXK (E) 2013-0069, respectively. After 1 week of normal feeding, high-fat and high-sugar diets and regular diets were prepared according to the formula. Twenty rats were fed with high-fat and high-sugar diet, high-fat and high-sugar diet feed formula: 10% lard, 20% sucrose, 2% cholesterol, 60% ordinary feed, 8% egg yolk powder(Jiangsu Xietong Pharmaceutical Bio-engineering Co., Ltd.) and 10 rats in normal control group were fed with normal diet. Twenty rats were fed a 100 g high-fat and high-sugar diet daily with free access to food and water. Corn oil 5 mL/kg was given by gavage on an empty stomach at 8:00 am every morning. At the end of the 12th week, 10 rats in the model group were given intraperitoneal injection of streptozotocin (STZ) (30 mg/kg) overnight fasting, and the other 10 rats in the model group were not given STZ injection. Ten control rats fed a regular diet were injected with the corresponding volume of citrate buffer. Blood samples were collected from the tail vein of 10 rats fed a high-fat and high-sugar diet and injected with STZ 3 days later, and random blood glucose was≥16.7 mmol/L. Blood glucose was monitored after 2 weeks, and fasting blood glucose was detected on the 2nd, 4th, 6th, 8th, 10th, 12th and 14th days. At the 14th week, 4 rats in each group were randomly selected to complete liver ultrasound, and fatty liver results were formed to determine the success of modeling. Finally, the rats were divided into 3 groups: (1) T2DM+MASLD group (n=10); (2) High-fat and high-glucose group (n=10); (3) normal control group (n=10); After 15 weeks, the rats were sacrificed under anesthesia, and the livers were collected, rinsed with normal saline at 4°C, weighed and placed in preservation solution and stored in a refrigerator at −80°C.

Real-time fluorescence quantitative PCR (Polymerase Chain Reaction): 0.15 g of rat liver tissue was taken from each group, and total RNA was extracted by the Trizol method. cDNA was reverse-transcribed into the corresponding cDNA in accordance with the reverse-transcription kit, and then subjected to real time Polymerase Chain Reaction (PCR)、RELB、S100A9 and SOCS1 primer sequences are shown in, Chain Reaction), PCR reaction total system 20μL, PCR amplification conditions: denaturation at 95°C for 10 min, annealing at 60°C for 1 min, extension at 95°C for 15 s, 40 cycles, each sample set up 3 replicate wells, internal reference GAPDH. The results of the experiments were analyzed by Bio-Rad Fluorescence Quantitative Analysis. The results of the experiments were read by Bio-Rad fluorescence quantitative analysis software, and the quantitative calculation of the mRNA expression levels of and in each group was expressed by 2−ΔΔCt (CT is the number of cycles). The main observation indexes were RELB, S100A9, SOCS1 mRNA expression levels in rat liver tissues in each group. 9

Kidney tissue specimens were fixed for 48h and paraffin-embedded. The paraffin-embedded tissue sections were serial sectioned at 3 μm, dewaxed and hydrated. 3% H 202 was incubated at room temperature for 10 min to remove endogenous peroxidase. After 10 min of citrate antigen repair, primary antibodies RELB (1:50), S100A9(1:300), SOCS1 (1:300), were added dropwise and incubated at 37°C for 25 min, and then freshly prepared DAB color development solution was added dropwise. The film was stained with hematoxylin and sealed with neutral gum. At the same time, PBS was used as negative control instead of primary antibodies. The immunohistochemical pathology pictures were analyzed using ImageScope software as follows: three different 400x fields of view were intercepted from each section, enter the analysis module of lmageScope software, and set all the dark brown color on the tissue sections as strong positive, brownish yellow as moderate positive, light yellow as weak positive, and blue nuclei as negative. Each tissue point was then identified and analyzed to find out the area (in pixels) of strong positive, moderate positive, weak positive and negative, the percentage of positivity and finally the H-score rating.

High-throughput sequencing expression profile data (GSE136103) were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/↗) (29). Twenty liver tissue samples were selected, including 9 MAFLD patients and 11 healthy individuals. Doublets in each sample were identified using DoubletFinder (30). Subsequently, doublets were filtered out, retaining cells meeting the following criteria: nFeature_RNA > 500 and < 3000; nCount_RNA < 10,000; percent_mito < 10. Cells from different samples were then merged using the merge function and analyzed using the Seurat package (31). Using the NormalizeData function, data were normalized for each library with the LogNormalize method and scale.factor of 10,000. The top 2,000 variable features were identified through the FindVariableFeatures function (32), and linear dimensionality reduction was performed using the RunPCA function with default parameters (npcs = 50). Batch effects were corrected using the RunHarmony function (33), followed by nonlinear dimensionality reduction using the RunUMAP function. The nearest neighbor graph was constructed using the FindNeighbors function (k.param = 20) with 30 principal components. Cell clustering was completed using the FindClusters function. Cluster marker genes were identified through the FindAllMarkers function with min.pct = 0.5 and logfc.threshold = 1. Single-cell samples were clustered and characterized using known lineage markers and manually annotated with reference to the human liver cell atlas on the CellMarker2 website (http://117.50.127.228/CellMarker/↗) (34). Differential expression analysis between cell clusters was performed using the FindMarkers function (Wilcoxon rank-sum test, p_val < 0.05 & abs(avg_log2FC) > 0.5).

To further explore intergroup heterogeneity across cell clusters, we analyzed proportional differences in cell counts between control and disease groups within distinct clusters, intergroup density variations, and distributional biases across conditions. Additionally, to test changes in cellular abundance at high resolution between groups/conditions, we employed the neighborhood abundance method in Milo (v2.2.0) (35) using the functions "buildGraph" and "makeNhoods" with parameters: "K=30", "d=50", "prop=0.2", "refined=TRUE", and "refinement_scheme=reduced_dim", indicating sampling using dimensionality-reduced data, with remaining parameters selected according to Milo's standard workflow recommendations. For pathway activity scoring of each cell cluster, we used the R package AUCell (v1.28.0): first, gene expression rankings per cell were computed via AUCell_buildRankings using the expression matrix with default parameters; then, custom-defined gene sets constructed from key genes identified in bulk-RNA analyses were applied to score individual cells; during this process, for each cell, the Area Under the Curve (AUC) value was calculated with AUCell_calcAUC based on gene expression rankings, where the AUC value represents the proportion of genes in the custom-defined set appearing among top-ranked genes per cell (36). Gene Set Variation Analysis implemented in the GSVA package (v2.0.0) was applied for gene set enrichment analysis using the "HALLMARK gene set" exported via GSEABase (v1.44.0). Intergroup differences in pathway activity scores per cell cluster between control and disease groups were computed using the LIMMA package (v3.62.2). Finally, based on intergroup heterogeneity patterns across cell clusters, we extracted Liver Sinusoidal Endothelial Cells (LSECs) and repeated the standard Seurat V5 workflow on these cells to identify and visualize endothelial subpopulations exhibiting significant distributional differences between groups.

Unsupervised pseudotemporal analysis was performed using the "Monocle" package (v2.34.0) (37). First, a cell dataset containing the expression matrix, phenotypic data, and feature data was constructed with the newCellDataSet function (expressionFamily = negbinomial.size). Next, size factor dispersion and gene expression across cells were corrected via estimateSizeFactors and estimateDispersions functions. Subsequently, dimensionality reduction was conducted using the DDRTree method (parameter max_components = 2), and cell ordering and visualization were performed with plot_cell_trajectory (38), while filtering highly variable genes correlated with pseudotime; expression changes of these genes along pseudotime were analyzed using the differentialGeneTest function. Finally, the filtered genes were clustered into distinct groups based on expression patterns and visualized using plot_pseudotime_heatmap. To identify genes bifurcating cells into different branches while exploring specific roles of key genes from bulk RNA analyses in single-cell pseudotime trajectories, we conducted "Branched Expression Analysis Modeling" (BEAM) analysis (39), further enhancing the filtering threshold; genes from BEAM analysis were visualized via plot_genes_branched_heatmap. Subsequently, KEGG functional enrichment analysis was performed for genes within each cluster using clusterProfiler (v4.14.4) and org.Hs.eg.db (v3.20.0) packages.

Initially, to further explore the LSEC clusters derived from prior research that may be closely associated with MAFLD progression, we processed and re-clustered LSECs using the standard Seurat V5 workflow. Then, by comparing these clusters through UMAP visualization of MAFLD and healthy control groups, the LSEC subpopulation exhibiting the most distinct cellular distribution was identified. Additionally, beyond LSECs, we identified cell clusters specifically distributing "bulk-RNA key genes" through scatter plots of gene distribution across cells; these cell clusters were simultaneously analyzed with LSECs for cell-cell communication analysis using the "CellChat" package (v2.1.2) (40). A random seed was set (seed=0528), and 8000 cells were randomly selected to create a CellChat object. Ligand-receptor interactions of "secretory signaling" were analyzed using "human" data from CellChatDB. The netVisual_circle function displayed the number and strength of inter-cell connections, while the netVisual_bubble function illustrated communications between different cell types. Subsequently, the netVisual_aggregate function demonstrated the communication network of a specific signaling pathway and calculated the contribution of various ligands. The mutual communication network between cells was visualized using the netVisual_individual function. Network centrality analysis was conducted using the netAnalysis_computeCentrality function (41) and represented as a heatmap.

Data analysis utilized SPSS 26.0 statistical software. Normally distributed measurement data were expressed as mean ± standard deviation and compared between two groups using independent samples t-test. Non-normally distributed measurement data were presented using the median and interquartile range, and comparisons between groups were conducted using the Mann-Whitney test. Statistical significance was set at P < 0.05.

The datasets GSE15653 and GSE89632 were downloaded from GEO data. The GSE89632 dataset was divided into GSE89632_SS group and GSE89632_NASH group according to SS and NASH samples. PCA principal component analysis was performed separately, and in the GSE15653 dataset, PC1 was 26.2% and PC2 was 8.7% (Figure 1A); in the GSE89632_SS dataset, PC1 was 30.6% and PC2 was 12.1% (Figure 1B); in the GSE89632_NASH dataset, PC1 was 30.3% and PC2 was 12.1% (Figure 1C). The volcano plot showed that 1686 differentially expressed genes were identified in the GSE15653 dataset with |log2(FC)|>0.5, p-value<0.05 as the screening threshold, of which 864 were up-regulated and 822 were down-regulated (Figure 1D). In the GSE89632_SS dataset, 2219 differentially expressed genes were identified, of which 1003 genes were up-regulated in expression and 1216 genes were down-regulated in expression, as in (Figure 1E). In the GSE89632_NASH dataset, 2248 differentially expressed genes were identified, of which 1062 genes were up-regulated in expression and 1186 genes were down-regulated in expression, as in (Figure 1F). The ring value heatmap and complex value heatmap showed the genes ranked top40 in differential folds in GSE15653 (Figure 1G), GSE89632_SS (Figure 1H) and GSE89632_NASH (Figure 1H), respectively.

Weighted gene co-expression network analysis (WGCNA) was performed on GSE89632SS, GSE89632NASH. Hierarchical clustering was performed based on similarity or correlation between samples and visualized as a hierarchical dendrogram (Figures 2A, E). The distance between samples indicated the degree of similarity or correlation between samples, and outlier samples were removed by pruning operation, and the outlier samples rejected were GSM2385767 and GSM2385782, and the most representative subset of samples was selected, and the clustering differences between different samples in the Normal group and the MAFLD group were visualized in a joint heat map (Figures 2C, G). The top 25% of genes with the largest fluctuations were selected for WGCNA analysis based on standard deviation ranking. The scale-free fit indices and average connectivity of various soft-threshold powers were evaluated on the basis of scale-free R2. Among them, GSE89632_SS selected soft threshold powers with β = 11 and scale-free R2 = 0.8 (Figure 2B), and GSE89632NASH selected soft threshold powers with β = 5 and scale-free R2 = 0.8 (Figure 2F). Gene clustering tree with module identification demonstrates the results of hierarchical clustering of genes in the two datasets and categorizes genes into different modules, further merging similar modules to reduce redundancy (Supplementary Figures 2A, D). The hierarchical clustering results between different modules were visualized using module clustering dendrogram (Supplementary Figures 2B, E), while the correlations between different modules were visualized with modal correlation heatmap (Supplementary Figures 2C, F). Finally, GSE89632SS identified 7 modules (Figure 2D), GSE89632_NASH identified 7 modules (Figure 2H). We performed KEGG enrichment analysis on genes within the "brown" modules of GSE89632_SS and GSE89632_NASH. The results demonstrated that in the SS group, pathways including TNF-JNK signaling pathway,IL-1/IL-1R-JNK signaling cascade, and TLR2/4-MAPK signaling axis collectively indicate the role of pro-inflammatory signaling in driving steatosis progression toward fibrosis. Concurrently, pathways such as Environmental factor-induced PI3K signaling pathway and Environmental factor-mediated RAS-ERK signaling cascade (activated by heavy metals) suggest synergistic interactions between environmental carcinogens and metabolic dysregulation. In the NASH group, the recurrent enrichment of TNF-JNK signaling pathway and Environmental factor-mediated RAS-ERK signaling cascade further validates conserved mechanisms underlying inflammation-carcinogenesis transition. The observed negative correlation may reflect compensatory dysregulation of protective mechanisms during disease progression, where downregulation of module genes attenuates their inhibitory effects on pro-inflammatory/carcinogenic pathways, thereby accelerating pathological advancement (Supplementary Files 1, 2). Based on this observation, we selected the "brown" modules from both SS and NASH WGCNA cohorts.

UpSet plots visualize the results of intersections of DEGs of GSE15653 taken with GSE89632_SS, GSE89632_NASH, SS WGCNA, NASH WGCNA, respectively. (Figure 3A). In order to better understand the interactions between genes, we used the STRING database and constructed the PPI network (Figures 3B, C), while the CytoHubba plugin was used to identify the top 15 hub genes in the GSE89632_SS group and the GSE89632_NASH group, respectively, based on the MCC, MNC, and Degree algorithms (Figures 3D–F, 3H–J), and then these three algorithms took the intersection to obtain 13 hub genes each (Figures 3G, K). Next, we investigated the chromosomal localization (Figures 3L, M) and interrelationships (Figures 3N, O) of the 13 hub genes.

We employed Lasso, Boruta, and RFE-SVM algorithms to identify hub genes associated with T2DM and MAFLD in the GSE89632_SS and GSE89632_NASH datasets. In GSE89632_SS, Lasso analysis revealed four significant variables: CLDN7, IRAK3, TFRC, and TNFRSF1A (Figures 4A, B). The Boruta algorithm identified ten significant variables, including CDKN1A, GJA1, IL1B, IRAK3, LBP, SOCS1, STAT3, TFRC, TLR2, and TNFRSF1A (Figure 4C). RFE_SVM highlighted eight features: IRAK3, MMP9, TNFRSF1A, TLR2, JUNB, TFRC, RELB, and CDKN1A (Figure 4D). In GSE89632_NASH, Lasso analysis identified five significant variables: ASPM, CCNF, CLEC4E, CX3CR1, and JUNB (Figures 4F, G). Furthermore, the Boruta algorithm filtered eleven significant variables, specifically ASPM, CCNF, CLEC4E, CX3CR1, GJA1, IL1B, IRAK3, JUNB, SOCS1, TLR2, and TNFRSF1A (Figure 4H). Two features were identified through RFE_SVM: JUNB and CX3CR1 (Figure 4I). Subsequently, we intersected the machine learning results from both datasets with autophagy-related genes to derive diagnostic biomarkers associated with autophagy in type 2 diabetes mellitus combined with NAFLD (Figures 4E, J).

We assessed the diagnostic performance of IRAK3, TFRC, TNFRSF1A, CX3CR1, and JUNB using the datasets GSE89632_SS, GSE89632_NASH, and GSE15653 for training, followed by plotting ROC curves. The results indicated that in GSE89632_SS, the AUC for IRAK3, TFRC, and TNFRSF1A exceeded 0.8 (Figure 4K). In GSE89632_NASH, the AUC for CX3CR1 and JUNB also surpassed 0.8 (Figure 4L), reflecting their high predictive accuracy. In GSE15653, the AUCs for all five genes were above 0.8 (Figure 4M). Subsequently, we conducted external validation using independent datasets GSE24807 and GSE23343. The results showed that in GSE24807, the AUC values of IRAK3, TNFRSF1A, CX3CR1, and JUNB were all greater than 0.6, while in GSE23343, the AUC values of IRAK3, TNFRSF1A, CX3CR1, and JUNB were all greater than 0.5 (Figures 4N, O). These results indicate that in the two independent external validation sets, these indicators also possess certain diagnostic value, which is consistent with our previous predictions. Ultimately, we identified IRAK3, TNFRSF1A, CX3CR1, and JUNB as diagnostic biomarkers for T2DM complicated with MAFLD associated with autophagy and endoplasmic reticulum stress.

During the GO and KEGG enrichment analysis, TNFRSF1A and JUNB were predominantly enriched in the TNF signaling pathway, while TNFRSF1A and TFRC showed significant enrichment in the HIF-1 signaling pathway. Furthermore, TNFRSF1A demonstrated substantial enrichment in the NF-kappa B signaling pathway, MAPK signaling pathway, Fluid shear stress and atherosclerosis, and Insulin resistance (Figure 5A, Figure 5C). CX3CR1 exhibited notable enrichment in the Cytokine-cytokine receptor interaction pathway, and TFRC displayed significant enrichment in the Hematopoietic cell lineage (Figures 5B, D).

In Enrichment analysis of GSEA, IRAK3、TNFRSF1A 和 JUNB commonly enriched in REACTOME_SIGNALING_BY_INTERLEUKINS(Supplementary Figure 3A, Supplementary Figure 3E, Supplementary Figure 3M).IRAK3 mainly enriched in BIOCARTA_IL1R_PATHWAY(Supplementary Figure 3B);REACTOME_TOLL_LIKE_RECEPTOR_TLR1_TLR2_CASCADE(Supplementary Figure 3C);REACTOME_TOLL_LIKE_RECEPTOR_CASCADES (Supplementary Figure 3D).TNFRSF1A mainly enriched in KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION (Supplementary Figure 3F);REACTOME_INTERLEUKIN_10_SIGNALING (Supplementary Figure 3G);KEGG_MAPK_SIGNALING_PATHWAY(Supplementary Figure 3H).TFRC mainly enriched in PID_HIF1_TFPATHWAY(Supplementary Figure 3I);KEGG_HEMATOPOIETIC_CELL_LINEAGE (Supplementary Figure 3J). JUNB mainly enriched in REACTOME_INTERLEUKIN_4_AND_INTERLEUKIN_13_SIGNALING(Supplementary Figure 3K);WP_TGFBETA_SIGNALING_PATHWAY(Supplementary Figure 3L); PID_IL6_7_PATHWAY(Supplementary Figure 3N);REACTOME_SIGNALING_BY_TGF_BETA_RECEPTOR_COMPLEX (Supplementary Figure 3O).

To assess immune cell infiltration and functional differences between the T2DM and control groups, we utilized the ssGSEA method to evaluate various immune cell subpopulations. The group comparison plot revealed increased levels of Activated dendritic cells, Monocytes, and Regulatory T cells in the T2DM group, with decreased levels of Memory B cells (Figure 6A). The correlation coefficient's absolute value indicates the strength of correlation: 0.3-0.5 denotes weak correlation, 0.5-0.8 signifies moderate correlation, and 0.8–1 represents strong correlation; P < 0.05 indicates statistical significance. In the T2DM group, IRAK3 exhibited a negative correlation with Natural killer cell infiltration level (R=-0.933) (Figure 6B, Supplementary Figure 4A). The correlations between TNFRSF1A and Effector memory CD8T cells, Gamma delta T cells, Macrophages, Natural killer cells, Type 1 T helper cells, and Regulatory T cells were positive (R=0.733, 0.700, 0.717, 0.717, 0.700, 0.733) (Figure 6C, Supplementary Figures 4B–G). CX3CR1 showed a negative correlation with Type 17 T helper cell infiltration (R=-0.983) (Figure 6D, Supplementary Figure 4H). JUNB displayed a positive correlation with Monocyte infiltration (R=0.700) (Figure 6E, Supplementary Figure 4I). Additionally, correlations existed between different types of immune cells (Supplementary Figure 4J).

The results of real-time fluorescence quantitative PCR analysis showed that there were significant differences in the expression levels of IRAK3, TNFR1, CX3CR1 and JUNB between the two groups (P < 0.05); among them, the mRNA expression levels of IRAK3, TNFR1, CX3CR1 and JUNB in the liver tissues of rats in the T2DM combined with MAFLD group were significantly higher than that in the control group (P > 0.05).), as shown in Figure 7F.

Immunohistochemical picture analysis showed that the expression levels of the four factors, CX3CR1, IRAK3, JUNB, and TNFR1, were significantly different between the T2DM & MAFLD groups and the Control group (P < 0.05); and the T2DM combined with MAFLD group were significantly higher than the Control group. For details, see Supplementary Table S3 and Figures 7A–E.

A comprehensive scRNA-seq analysis was conducted using publicly available data from healthy and MAFLD liver samples. The clustering results, illustrated through dendrograms(Supplementary Figure 5A) and marker genes bubble plots (Figure 8B), culminated in the identification of 13 major cell clusters at a resolution of 0.6. Distinct cell types were pinpointed, including Cycling(cycling cells) (STMN1,MKI67), B_cell (CD79A,CD79B,MS4A1), Plasma(Plasma cell) (CD79A,IGHA1), pDC (plasmacytoid dendritic cell) (LILRA4,CLEC4C), Mes(mesenchymal cell) (PDGFRB,ACTA2,COL1A1,COL1A2,COL3A1,DCN), Cho (Cholangiocyte) (EPCAM,KRT19,FXYD2), T_cell (CD3D,CD3E,CD3G,CD8A), ILC (innate lymphoid cell) (KLRF1,KLRC1,GZMB,NKG7), Kupffer (CD163,MARCO), Neu(Neutrophils) (S100A8,S100A9), cDC(Conventional Dendritic Cell) (CLEC10A,CD86), LSEC(Liver Sinusoidal Endothelial Cel) (CLEC14A,CD34,VWF,STAB2,CLDN5) and LEC(Lymphatic Endothelial Cells) (PROX1,TFF3,CCL21). The findings were further visualized using UMAP (Figure 8A).

We analyzed the proportional differences in cell counts between control and disease groups across distinct cell clusters, with LSEC showing a significantly higher proportion in the MAFLD group (Figure 8C), along with proportional variations of cell clusters among different samples (Figure 8D). We performed density analysis of inter-group cell counts at a two-dimensional level, where relative brightness indicates higher density of the corresponding cell cluster in its respective group (Figures 9A, B). Subsequently, we evaluated the distributional preference between MAFLD and Control groups using the Ro/e index (Figure 9C), as well as distributional preferences across individual samples (Figure 9D), revealing LSEC exhibiting stronger enrichment tendencies in the disease group. Further, in the Milo differential abundance testing, we identified significant disease-state-driven changes occurring specifically in LSEC, with a negative Log Fold Change in relative abundance (Figures 9E, F). Despite LSEC's preferential enrichment in the MAFLD group (Figure 9C), the actual abundance of LSEC cells was significantly reduced (Figure 9F), suggesting potential widespread LSEC activation or proliferation in MAFLD state alongside intrinsic heterogeneity within LSEC populations, where certain subpopulations may undergo suppression or reduction. Subsequently, we explored cellular heterogeneity through differential gene expression and enrichment scoring perspectives to elucidate connections between key genes and various cell clusters. Using thresholds of |avg_log2FC|>1 and p_val<0.05, the FindMarkers function identified differentially expressed genes suitable for a multi-group volcano plot (Figure 10A). The potential marker genes for diabetes-combined MAFLD, derived from bulk RNA analysis, displayed significant differences during single-cell differential analysis. Specifically, JUNB was upregulated in pDC, Neu, and kupffer, but downregulated in LEC; CX3CR1 was upregulated in Neu, LSEC, kupffer, ILC, and cDC; IRAK3 exhibited downregulation in Neu and LSEC. CX3CR1 (Figures 10C, D) and IRAK3 (Figures 10E, F) also showed distinct distributional variations across cell types, predominantly localized within five cell populations: LSEC, ILC, Neu, cDC, and kupffer. Additionally, using key genes including JUNB, CX3CR1, IRAK3, and TNFRSF1A to construct a custom gene set, we scored activity levels across cell clusters on this custom dataset. Results revealed relatively higher activity of the key gene-defined custom set within these five cell types: Neu, LSEC, kupffer, ILC, and cDC (Figure 10B). Finally, we focused our investigation on LSEC and the other four cell types (ILC, Neu, cDC, kupffer) for subsequent research. Through gene set variation analysis in LSEC, we enriched pathways including "TGF BETA SIGNALING", "WNT BETA CATENIN SIGNALING", and "HEDGEHOG SIGNALING" (Figure 11A). The activity of these pathways indicates that in MAFLD, LSECs likely reside in a highly dynamic microenvironment involved in regulating cell proliferation, differentiation, immunity, and cellular interactions. When comparing enrichment differences in LSEC between disease and MAFLD groups, we observed significantly higher enrichment of "HEDGEHOG SIGNALING" in the MAFLD group. Additionally, pathways including "MYC TARGETS V2", "PI3K-AKT-MTOR SIGNALING", "IL6-JAK-STAT3 SIGNALING", "APOPTOSIS" and "PROTEIN SECRETION"—significantly associated with biological processes like cell proliferation, survival, metabolic regulation, immune response, inflammatory response, and programmed cell death—exhibited substantial inter-group enrichment differences (Figure 11B). These results demonstrate that LSECs likely play crucial roles in MAFLD pathogenesis, warranting further investigation into their specific mechanisms during MAFLD development.

Re-clustering of LSEC through secondary dimensionality reduction yielded 13 subpopulations, among which clusters 2, 6, 7, 8, and 12 showed significant reduction compared to the control group. This aligns with prior Milo analysis conclusions, indicating intrinsic heterogeneity within LSEC populations where specific subpopulations may be suppressed or diminished in the disease state. Conversely, clusters 0, 1, 5, 9, and 10 exhibited marked expansion. We defined the expanded clusters 0, 1, 5, 9, 10 as "LSEC_inc" and the reduced clusters 2, 6, 7, 8, 12 as "LSEC_red" (Figure 12A). We conducted a pseudo-temporal trajectory analysis on LSEC, dividing the entire trajectory into five stages (Supplementary Figure 5B), showing the distribution of different LSEC subgroups in the trajectory (Supplementary Figure 5C) and the direction of differentiation and development (Supplementary Figure 5D). The cell density plot along the timeline further illustrates the distribution and dynamic changes of LSEC during the pseudo-temporal process (Supplementary Figure 5E). There is a significant difference in the temporal trajectory distribution of LSEC between the Control group and the MAFLD group (Supplementary Figure 5F). Subsequently, we analyzed the expression patterns of highly variable genes in LSEC during the fitting process. These highly variable genes were clustered into four different clusters. Among our key genes of interest, JUNB—identified as a highly variable gene along the pseudotime trajectory—was grouped into the fourth cluster (Figure 12B). To further explore the role of JUNB in pseudotime progression, through BEAM (Branched Expression Analysis Modeling) analysis, we identified the potential regulatory genes that regulate the branches of LSEC in the pseudo-temporal process, and clustered these genes into four different clusters according to different expression patterns. In BEAM analysis, JUNB was also clustered as a highly variable gene in the first cluster (Figure 12C). These results suggest that JUNB may participate in critical regulatory processes during pseudotime. To investigate further, we performed KEGG enrichment analysis on genes in the first cluster (containing JUNB) from BEAM analysis. Consequently, we enriched pathways potentially relevant to LSEC, including "Fluid shear stress and atherosclerosis" and "Leukocyte transendothelial migration" signaling pathways, indicating that LSEC alterations during MAFLD progression may relate to fluid shear stress and leukocyte transendothelial migration (Supplementary Figure 5G, Supplementary Table S4).

We employed the CellChat package to analyze cell-cell communication networks among five key cell types in MAFLD: LSEC, ILC, Neu, cDC, and kupffer. "LSEC_red" and "LSEC_inc" were defined as cell types at an equivalent hierarchical level to ILC, Neu, cDC, and kupffer. Circular plots depict interaction quantities (Supplementary Figure 6A) and interaction strengths (Supplementary Figure 6B) among cellular populations. When "LSEC_inc" acts as a signal receiver, it communicates with multiple cell types through ligand-receptor pairs including CCL5-ACKR1 and PPIA-BSG (Figure 13A); when serving as a signal sender, it interacts via MIF-(CD74+CD44), MIF-(CD74+CXCR4), and PPIA-BSG (Figure 13B). For "LSEC_red" as a receiver, communications occur through GZMA-F2R and PPIA-BSG pairs (Figure 13C); as a sender, interactions involve MIF-(CD74+CD44), MIF-(CD74+CXCR4), and PPIA-BSG pairs (Figure 13D). The dominant communicating ligand-receptor pairs remained consistently identified, with corresponding signaling pathways classified as: MIF-(CD74+CD44) and MIF-(CD74+CXCR4) belong to the MIF signaling pathway; CCL5-ACKR1 belongs to the CCL signaling pathway; GZMA-F2R belongs to the PARs signaling pathway; PPIA-BSG belongs to the CypA signaling pathway (Supplementary Tables S5, S6). Consequently, this study focuses on the aforementioned four signaling pathways.

Within the MIF signaling pathway, both "LSEC_inc" and "LSEC_red" act exclusively as signal senders when communicating with other cells, displaying no self-communication. In hierarchical diagrams corresponding to MIF-(CD74+CD44) and MIF-(CD74+CXCR4), extensive communications exist between LSEC and ILC, Neu, cDC, and kupffer (Figures 13E, 14A). In the CypA pathway, both "LSEC_inc" and "LSEC_red" communicate with kupffer cells as signal senders (Figure 14B). Conversely, within the PARs signaling pathway, nearly exclusive communication occurs where "LSEC_red" acts as a signal receiver interacting with ILC (Figure 14C). The CCL5 signaling pathway exhibits an inverse pattern, with "LSEC_inc" primarily functioning as a signal receiver communicating with ILC (Figure 15A). Subsequently, we ranked ligand-receptor pair contributions across these four pathways (Supplementary Figures 6D–G) and consolidated the results. Compared to other ligand-receptor pairs, two pairs from the MIF pathway—"MIF-(CD74+CXCR4)" and "MIF-(CD74+CD44)"—along with the PPIA-BSG pair from the CypA pathway demonstrated significantly higher pathway contributions (Supplementary Figure 6C). This highlights the critical importance of these two signaling pathways and the three key ligand-receptor pairs in mediating communications between LSEC and ILC, Neu, cDC, and kupffer.

Finally, by computing network centrality scores, we identified primary signaling roles of distinct cells within the specified pathways mentioned above (Figures 15B–E, Supplementary Figure 7). Additionally, analysis of signaling roles across all aggregated cell-cell communication networks revealed comprehensively assessed communication identities for different cell types. Both MIF and CypA signaling pathways exhibited the highest activity among all pathways, which aligned with our pathway contribution ranking (Figure 15F).