The macrophage-associated microRNA-4715-3p / Gasdermin D axis potentially indicates fibrosis progression in nonalcoholic fatty liver disease: evidence from transcriptome and biological data

ABSTRACT Nonalcoholic fatty liver disease (NAFLD) is highly possible to progress to cirrhosis, malignancy, and liver failure through fibrogenesis. The enormous potential of pathogenetic and therapeutic targets in NAFLD has been revealed. This study aimed to explore novel factors potentially indicating or mediating NAFLD progression. Multiple bulk and single-cell RNA sequencing datasets were used, in which landscapes of cell populations were clarified to characterize immune cell infiltration. Significantly high infiltration of macrophages (MPs) was discovered during NAFLD progression. Samples in bulk NASH datasets were regrouped by MP level. Highly differentially expressed genes (DEGs) were identified in the Ctrl vs. NASH comparison, low MP vs. high MP comparison, and the weighted gene co-expression network analysis (WGCNA) clusters. Eight hub genes were identified as promising targets by protein–protein interaction analysis and validated in fibrosis progression, microRNA (miR)–protein interactions were predicted, and the hub genes were verified in a free fatty acid (FFA)-induced macrophage injury model. The results showed that Gasdermin D (GSDMD) was upregulated with fibrosis progression in NAFLD and was associated with macrophage infiltration. In addition, a potential regulator (miR-4715-3p) was correlated with GSDMD. The miR-4715-3p/GSDMD axis potentially modulates macrophage-associated immunity and indicates fibrosis progression in NAFLD.


Introduction
onalcoholic fatty liver disease (NAFLD) is one of the most common chronic liver diseases worldwide. Lipid deposition and steatosis are independent causes of NAFLD, independent of drug-related, alcohol-related, or genetic risks. The pathogenesis and progression of NAFLD are fundamentally mediated by the complex interactions among steatosis, fibrosis, and inflammation [1]. Although early symptoms rarely appear in nonalcoholic steatohepatitis (NASH) patients, chronic liver cirrhosis eventually leads to liver failure or hepatocellular carcinoma (HCC). Notably, multiple types of immune modulation occur during NAFLD progression. As executors and directors in the immune microenvironment, immune cells, especially macrophages and T cells, are requisite factors in the pathogenesis, cirrhogenesis, metabolic compensation, and progression of NAFLD [2,3].
In almost all chronic liver diseases, fibrosis occurs and progresses with inflammatory reactions, which are involved in liver macrophage modulation. Systematically, the liver macrophage population generally consists of Kupffer cells (KCs) and monocyte-derived macrophages (MoMFs). Liver macrophages play an essential role in fibrogenesis through interactions with hepatic stellate cells (HSCs) [4]. As a feedback regulatory mechanism, chemokines and cytokines secreted by HSCs enhance macrophage infiltration and expansion, which promote the fibrotic phenotypes and survival of HSCs. Passively activated M2 macrophages are associated with liver injury in NAFLD and mediate fibrotic responses conducive to liver remodeling and regeneration by secreting transforming growth factor-β (TGF-β) and platelet-derived growth factor (PDGF) [5]. Although crucial roles of macrophages in NAFLD and other inflammatory liver diseases have gradually been revealed, the exact mechanisms of macrophage involvement in the pathogenesis and progression of these diseases remain unclear. In recent decades, transcriptome analysis was rapidly increased and has provided logical and scientific guidance to traditional biological studies. Fundamentally, computational analysis based on high-throughput datasets can indicate both target genes and biological processes regarding specific research goals [6,7]. Therefore, in NAFLD with a chronic fibrosis progression, bioinformatics analysis would bring out unique advantages to understand disease characteristics and therapeutic targets.
The hypothesis of this study is that miR-4715-3p/GSDMD axis associating macrophage infiltration potentially indicates NAFLD progression. Based on the analysis of high-throughput bulk/ single-cell RNA sequencing datasets, this study aims to reveal immune cell landscapes in NAFLD livers during fibrosis progression. Thereby, through further MP-associated analyses, the miR-4715-3p/GSDMD axis could be determined as a potential mediator in liver macrophages and cirrhosis progression.

Data resources and DEG analysis
In this study, expression profile datasets from bulk RNA-sequencing and single-cell RNA (scRNA)sequencing analyses were downloaded from the Gene Expression Omnibus (GEO) database (http:// www.ncbi.nlm.nih.gov/geo). Basic information is listed in Table 1. Probes from the raw data were matched with the official gene symbols with the DAVID online tool (https://david.ncifcrf.gov/) [8] and noncoding RNAs were excluded. Differentially expressed gene (DEG)-related analyses were performed using R Studio (limma package) [9]. Fold change (FC) >2 and p < 0.01 were used as the statistical criteria to identify significant DEGs. In addition, overlapping gene clusters were identified by the generation of Venn diagrams (http://bioinfor matics.psb.ugent.be/webtools/Venn/) [10].

Immune cell infiltration analysis and dataset regrouping
To investigate immune cell infiltration, raw data from bulk RNA-seq datasets were included and evaluated by multiple computational methods. The xCell package (https://xcell.ucsf.edu/) [11] in R Studio was applied, and the results were compared and validated with the CIBERSORT (https://cibersort.stanford.edu) [12] and EPIC (https://github.com/GfellerLab/ EPIC) [13] methods. Then, raw data from 74 NASH liver samples in the GSE164760 dataset were classified into the high MP and low MP groups based on the cell-type enrichment scores from Xcell. Thereby, DEG-related analyses were performed using R Studio (limma package).

Weighted gene co-expression network analysis (WGCNA)
To outline the gene expression patterns in multiple samples, the WGCNA package (http://www. g e n e t i c s . u c l a . e d u / l a b s / h o r v a t h / CoexpressionNetwork/Rpackages/WGCNA), which is a comprehensive method used to cluster co-expressed genes and compare the associations between modules and specific traits, was used with R Studio [14]. The expression matrix of NASH samples in the GSE164760 dataset regrouped by MP level was analyzed by WGCNA. The correlation analysis implemented in the WGCNA package is based on the Pearson method, and specific parameters (min Module Size = 30, reassign Threshold = 0, and merge Cut Height = 0.25) were used to run the program. The clusters of co-expressed genes were then displayed as modules labeled in multiple colors, with correlation coefficient values and significance (p values).

Cell culture and treatment
The human monocyte line THP-1 was purchased from the National Collection of Authenticated Cell Cultures (Shanghai, China) and cultured in RPMI 1640 medium supplied with 10% FBS (Gibco, Australia) in an incubator with 5% CO 2 at 37°C. Macrophage differentiation of THP-1 monocytes was induced with 100 nmol/L phorbol 12-myristate-13acetate (PMA) (Sigma-Aldrich, USA) for 48 h. Palmitic acid (PA) and oleic acid (OA) (both from Sigma-Aldrich, USA) were dissolved at a 2:1 ratio (final concentration of 30 µM) in RPMI 1640 medium supplemented with 0.1 M NaOH and 0.1% BSA (both from Sigma-Aldrich, USA) [25].

RNA quantification
Total RNA was extracted from cultured cells using TRIzol® reagent (Invitrogen; Thermo Fisher Scientific, USA). Reverse transcription was conducted with a PrimeScript™ RT Reagent Kit (Takara Biotechnology Co., Ltd.), and real-time quantification was subsequently performed using a SYBR® Premix Ex Taq Kit (Takara Biotechnology Co., Ltd.). The relative RNA expression levels were calculated using the 2 -ΔΔCt method [26]. The primer pairs used in this study are listed in Supplementary Table 1. SPSS 26.0 (SPSS Statistics, USA) was used for general statistical analysis. GraphPad Prism 9.0 software (GraphPad Software, USA) and R Studio were used to generate plots. Significant differences were identified using one-way ANOVA with the Bonferroni post hoc test. Pearson correlation analysis was performed to calculate correlation coefficients. Data are presented as the mean ± SEM values. p < 0.05 indicated a statistically significant difference (labeled '*').

Results
We hypothesized that miR-4715-3p/GSDMD axis may play as a mediator for macrophage infiltration and NAFLD progression. Multiple RNA-seq datasets were included in this study, and therefore, immune infiltration and hub genes were determined. GSDMD was identified as the key factor both associated with high macrophage infiltration and NAFLD progression. Predicted as a regulator to GSDMD, miR-4715-3p was validated in FFA-  induced macrophage injury, which was significantly correlated with GSDMD.

Identification of DEGs and immune cell landscape in NASH livers
The overall flow chart of the technical paths followed in this study is shown in Figure 1. In 74 NASH and 6 healthy control liver samples, DEGs were identified as genes with |log2 (fold change)| > 1. 5  KIAA1147 (DENND11), and TIPLR. The top 10 downregulated genes were PDIA4, SERPING1, CDK2AP2, C19orf60, NT5C, BLVRB, ASPA, MRPL41, GALR3, and ADRA1B. GSEA was then used to perform functional enrichment analysis on the upregulated gene cluster and visualize the results. Among the leading enrichment patterns, biological processes associated with immune cells (Th cells and macrophages) were revealed to be statistically significant (Figure 2c).
Accordingly, Xcell, CIBERSORT, and EPIC were applied to investigate and compare immune cell infiltration in NASH. Excluding nonimmune cells, the main myeloid cells and lymphocytes were evaluated and displayed. As shown in Figure 3a, the infiltration levels of B cells, Th2 cells, basophils, MPs, and M2-MPs were significantly increased in the NASH group compared with the control group. In contrast, the infiltration level of Th1 cells was decreased. As shown in Figure 3b, the infiltration levels of CD4+ T cells, T helper cells, gamma delta T cells, M0-MPs, and M2-MPs were significantly increased in the NASH group compared with the control group. In contrast, the infiltration level of plasma cells was decreased. As shown in Figure 3c, the infiltration levels of fibroblasts and MPs were increased. These findings conclusively show that MP infiltration is commonly elevated in NASH livers compared with healthy livers. To further explore the correlations between MP infiltration and cirrhosis progression in NAFLD, MP clusters (total MPs, M1-MPs, and M2-MPs) were evaluated with xCell in multiple RNA-seq datasets (GSE89632, GSE49541, and GSE139602). Increasing trends in infiltration were found for both steatosis (Figure 3d) and cirrhosis progression (Figure 3e and f).

The identification of macrophageassociated hub genes in NAFLD progression
Based on the xCell scores of MPs, 74 NASH liver samples from GSE164760 were divided into the high-MP (n = 37) and low-MP (n = 37) groups ( Figure 4a). The DEGs were then sorted and are shown in Figure 4b. In addition, WGCNA was conducted to generate gene coexpression modules, from which the module-trait relationships ( Figure 4c) and cluster dendrograms (Figure 4d) were generated. Considering both the correlation coefficients and p values, the MEblue module, with 3166 genes, was identified. After overlapping three gene clusters, 171 genes were identified as the most significantly upregulated macrophageassociated genes in NASH livers (Figure 4e). Then, GO enrichment analysis was carried out with these 171 genes, and the results are shown in Figure 4f; the terms NLRP3 inflammasome, interleukin-1 beta (IL-1β), and tumor necrosis factor (TNF) were found to be significantly enriched. A PPI network was constructed to visualize functional interactions. Hub clusters and seed genes were also identified ( Figure 4g).
To assess clinical correlations, the expression levels of eight hub/seed genes (GSDMD, REEP5, CHD9, TNRC6A, UBE3A, PGAM1, PAPOLA, and KPNA1) were verified on bulk RNA-seq datasets (GSE49541, GSE164760, and GSE139602) and single-cell RNA-seq datasets (GSE123661, GSE136103, and GSE98782). All eight genes were upregulated in NASH livers at advanced fibrosis stages compared with mild fibrosis stages (Figure 5a). In contrast to NASH livers, GSDMD and TNRC6A were upregulated in cirrhotic livers (Figure 5b). Among all stages of progression, a continuously increasing trend in GSDMD and TNRC6A expression was identified (Figure). At the single-liver MP level, GSDMD was significantly upregulated in cirrhotic KCs (Figure 5d). The distribution of cell populations in NAFLD-related cirrhosis is shown in Figure 5e. After single-cell analysis, GSDMD was found to be upregulated in the liver MP population (Figure 5f). Moreover, among macrophages in cirrhotic livers, monocytederived macrophages (MoMFs) expressed more GSDMD than KCs (Figure 5g). These results indicate that GSDMD may not only play crucial roles as a biomarker for NAFLD progression but also function through liver macrophages.

Correlative upregulation of miR-4715-3p and GSDMD in FFA-induced macrophage injury
By searching an online database, multiple miRNA clusters highly correlated with GSDMD were identified. Then, hsa-miR-4715-3p was found in all miRNA clusters (Figure 6a). To mimic immune modulation in NAFLD, an in vitro model of FFAinduced macrophage injury was applied. Then, THP-1 cells were treated sequentially with PMA to induce differentiation into macrophages and with FFA solution for 1 week (the cytotoxicity of FFA in THP-1-derived macrophages is shown in Supplementary Figure 1). MoMFs were harvested at sequential time points (on the first, second, third, and seventh days) for RNA quantification. The expression levels of hsa-miR-4715-3p ( Figure 6b) and the eight hub genes (Figure 6c) are shown at different time points. hsa-miR-4715-3p and GSDMD were generally upregulated after FAA treatment, and their expression peaked on the third day. In contrast, only REEP5 was upregulated at all time points. Therefore, Pearson correlation analysis was conducted on hsa-miR-4715-3p and the hub genes (Figure 6d), and GSDMD expression was found to be significantly positively correlated with hsa-miR-4715-3p expression (R = 0.96, p = 0.04). The results indicate that coexpression of miR-4715-3p and GSDMD is significantly induced in response to FFA stimulation in macrophages.

Discussion
NASH is considered an inflammatory subtype of NAFLD with steatosis and evidence of hepatocyte injury and multiple immune cell interactions [27]. Furthermore, fibrosis consistently occurs with hepatic cell injury and gradually progresses to cirrhosis, which results in hepatocarcinoma and liver failure [28]. Generally, steatosis and liver function can be measured by serological and imaging approaches. However, liver fibrosis could historically be assessed only by biopsy, which is sensitive and accurate but results in unavoidable injury and pain in patients [29]. Currently, the most effective noninvasive methods are ultrasonic transient elastography (for example, fiber scanning) and magnetic resonance elastography. Moreover, noninvasive estimation systems can evaluate the fibrosis stage in patients without biopsy, with the NAFLD fibrosis score and fibrosis-4 index (FIB-4) commonly used. The complex pathophysiological mechanism of NASH progression remains unrevealed [30,31]. Currently, multiple therapies aimed at several targets, such as changes in the microbiome and intestinal permeability, oxidative stress, insulin resistance, apoptosis, lipotoxicity, inflammation, bile acid metabolism, and fibrogenesis. Given that the therapeutic effects of traditional drugs on NASH progression are poor, the discovery of disease targets and research on targeted therapies are of great importance [32]. Generally, liver macrophages (including KCs and MoMFs) play pivotal roles in the progression and resolution of fibrosis. With the development of fibrosis, injuryinduced inflammation further enhances macrophage recruitment and HSC activation [33]. In addition, a variety of matrix metalloproteinases (MMP-9, MMP-12, and MMP-13) secreted by liver macrophages participate in matrix degradation, thus contributing to the alleviation of liver injury and fibrosis [34,35].
The initial goal of this study was to characterize immune cell infiltration and identify progressionrelated mediators in NAFLD. Bulk and single-cell RNA sequencing datasets of NASH and fibrosis progression were introduced as basic data sources. Then, multiple bioinformatic methods and databases were used to reveal the immune cell landscape and important factors. Of note, a novel method of dataset regrouping was developed in this study based on xCell and WGCNA to identify macrophage-associated DEGs. Subsequently, hub genes were identified from the PPI network. Further screening and validation were carried out with single-cell RNA-seq datasets and a model of FFA-induced macrophage injury. As a result, GSDMD and hsa-miR-4715-3p were found to be co-upregulated, indicating NAFLD progression.
GSDMD is a key pyroptotic substrate of inflammation-triggered caspase modules. GSDMD initiates pyroptosis directly via caspase-1-induced cleavage. Moreover, GSDMD-related pyroptosis mainly occurs in macrophages and is directly mediated by inflammasomes [36,37]. In addition to pyroptotic factors (caspase-1/-4 and IL-1β), NLRP3 and NLRC4, which act as classical regulators of inflammasomes [38], were also included in the final gene cluster. Accordingly, inflammasomeinduced pyroptosis is considered a crucial macrophage-fate-accompanying NAFLD progression. In addition, miRNA-protein interaction analysis revealed that hsa-miR-4715-3p might act as a potential upstream regulator of GSDMD. Previous studies have revealed multiple roles and function of several microRNAs in NAFLD and liver fibrosis [39][40][41][42]. According to a literature review, hsa-miR-4715-3p has been reported in only a few cancer studies [43,44]. In conclusion, the miR-4715-3p/GSDMD axis may be a novel signaling pathway that mediates not only macrophage functions but also NAFLD progression. Although this study provided robust evidence at the transcriptome level, further investigations are required to validate the functions of the GSDMD protein as an executor of inflammasome-mediated pyroptosis.
However, this study is mostly based on RNA-seq datasets and computational analysis, leading to inevitable limitations. Baseline levels and batch effects of different clinical cohorts could not be entirely adjusted, referring to statistical biases. Moreover, molecular validations were included but not emphasized in this study. Hence, to determine the exact relationship between GSDMD and hsa-miR-4715-3p, more molecular studies are required.

Conclusions
During fibrosis progression in NAFLD, the macrophage infiltration and highly associated gene -GSDMD -are elevated. The hsa-miR-4715-3p is upregulated in FFA-induced macrophage injury and parallelly altered by GSDMD expression. Thus, the macrophage-associated miR-4715-3p/GSDMD axis potentially indicates fibrosis progression in NAFLD.