Signature of gene expression profile of liver sinusoidal endothelial cells in nonalcoholic steatohepatitis

Background: There has been emerging evidence that liver sinusoidal endothelial cells (LSECs) play an important role in the pathogenesis of nonalcoholic steatohepatitis (NASH). This study aims to figure out the signature of the gene expression profile of LSECs in NASH and to explore potential biomarkers related to damaged LSECs in NASH.

Methods and materials: Animal experiments were performed to demonstrate the significant structural damage of LSECs in the NASH model. To further understand the functional changes of these damaged LSECs in NASH, we used the public GEO database that contained microarray data for the gene expression of LSECs in NASH and normal mouse liver. Differentially expressed genes (DEGs) were analyzed, and further Gene Ontology (GO) enrichment analysis was performed to understand the functional changes. The hub genes were then identified and validated via external GEO databases.

Results: There was significant structural damage to LSECs in the NASH model, accompanied by remarkable functional changes of LSECs with 174 DEGs (156 upregulated and 18 downregulated genes). The functions of these DEGs were mainly enriched in the inflammatory reactions and immune responses. Nine specifically expressed hub genes were identified. Among them, CCL4 and ITGAX showed the most significant correlation with NASH, with AUROC of 0.77 and 0.86, respectively. The protein–protein interaction network, mRNA–miRNA interaction network, and ceRNA network were further predicted.

Conclusion: LSECs show significant structural damage and functional changes in NASH. The LSEC-related DEGs, such as CCL4 and ITGAX, might be promising biomarkers as well as potential treatment targets for NASH.

Introduction

Nonalcoholic steatohepatitis (NASH) has become a leading cause of cirrhosis and hepatocellular carcinoma over recent years with increased overall mortality (Dufour et al., 2020). The ballooning degeneration of hepatocytes and lobular inflammation are the histological characteristics of disease progression from nonalcoholic fatty liver disease (NAFLD) to NASH (Kleiner et al., 2019). Early diagnosis of NASH is of critical importance to reduce the economic burden for patients with NASH. However, the liver pathological biopsy, the gold standard for diagnosing NASH, is still difficult to obtain for most patients.

The liver contains a large number of hepatocytes and multiple types of non-parenchymal cells, including liver sinusoidal endothelial cells (LSECs), Kupffer cells, natural killer (NK) cells, B lymphocytes, T lymphocytes, and hepatic stellate cells. Recent reports have brought the importance of LSECs to our attention as the role of the hepatic barrier against invasive factors (Wang and Liu, 2021a; Wang et al., 2021). However, under pathological conditions, LSECs can shift to a pro-inflammatory phenotype and thus contributing to NASH progression (Hammoutene and Rautou, 2019). LSECs are not considered innate immune cells, but they might show immunomodulatory functions under stressed conditions (Cai et al., 2019). For instance, in many types of liver diseases, LSECs can mediate the adhesion of circulating leukocytes (Carambia et al., 2014; Wang and Liu, 2021b). The damage of LSECs has been considered important for the activation of immune reactions in the progression of NASH. Injured LSECs show both structural and functional changes. For instance, in congestive hepatopathy, the mechanical stretch of LSECs could increase the expression of CXCL1 via integrin-dependent activation of transcription factors (Hilscher et al., 2019). LSECs can also release mechanosensitive angiocrine signals that further promote portal hypertension by recruiting sinusoidal neutrophils and promoting the formation of neutrophil extracellular traps and microthrombi (Hilscher et al., 2019). Under the stimulation of lipopolysaccharide (LPS), injured LSECs can express significantly increased chemokines and inflammatory factors and further mediate the recruitment of hepatic neutrophils to accelerate liver injury (Ng-Nguyen et al., 2021). The exact mechanisms of injured LSECs in inducing hepatic inflammation and liver diseases remain to be further investigated.

Currently, transcriptomic and microarray analysis has been widely applied to explore new biomarkers in a variety of diseases including cancers (Krishnan et al., 2021), metabolic diseases (Fujimoto et al., 2020), and auto-immune diseases (Cheng et al., 2021). In addition, verification of the competitive endogenous RNA (ceRNA) networks can present a deeper mechanism in the transcriptional regulatory network (Demircioglu et al., 2019). While a remarkable change in transcription is recognized as a defining feature of NASH, the potential key genes and pathway networks of LSECs are worth thorough analysis. Thus, in this study, we aim to comprehensively analyze the signature of the gene expression profile of LSECs in NASH.

Our study suggests great changes in the gene expression profile of LSECs in NASH and also predicts new biomarkers and pathways to further explore the interplay of regulatory mechanisms.

Methods

Animal experiments

Animal experiments were performed using 250–300 g male Sprague–Dawley rats, which were obtained from Beijing Vital River Laboratory Animal Technology Co., Ltd. The rats were fed with the methionine–choline-deficient diet (MCD) (obtained from Jiangsu Medicine Co., Ltd) for 6 weeks to induce NASH, while a normal chow diet was used for the control group. Rats were kept in cages in a temperature-controlled room (21 ± 2°C) with free access to food and water. Animal experiments were approved and supervised by the Ethics Committee of Peking University People’s Hospital (Approval no. 2020PHE074).

Measurement of alanine aminotransferase (ALT) and aspartate aminotransferase (AST)

At the appropriate time point of disease modeling, rats were euthanized by CO₂ asphyxiation, and blood samples were collected. After centrifugation of blood samples for 10 min at 3,200 rpm, serum was collected. ALT and AST were tested to evaluate liver injury by the assay reagent kits (Nanjing Jiancheng Bioengineering Institute, China) according to the manufacturer's instructions.

H and E staining

Liver tissue was fixed in paraformaldehyde and embedded in paraffin. The liver sample was cut into 4 μm thick sections and then baked at 60°C for 4 h and embedded in paraffin. After paraffin was removed by a graded ethanol series and xylene, sections were stained with hematoxylin and eosin.

Scanning electron microscopy (SEM) and transmission electron microscopy (TEM)

The perfusion fixation of rat liver tissue for light and electron microscopy was performed as previously described (Wisse et al., 1984). Briefly, under isoflurane anesthesia, the liver of SD rats was perfused to scour blood, followed by glutaraldehyde fixative buffer via the portal vein. The perfusion was stopped until the surface of the liver changed from red and soft to pale and tough. Liver samples were cut into 2 × 2 × 2 mm SEM (JSM-7900F, JEOL) and appropriate slices for TEM.

Microarray data acquisition

The public GEO database (https://www.ncbi.nlm.nih.gov/gds/) was searched to obtain microarray data for LSECs in NASH and normal mouse liver. The screening criteria included (1) expression profiling that was tested by array (2); LSECs were isolated from NASH and normal samples from liver biopsies (3); datasets contained no less than five samples (4); datasets contained full-scale information about the samples (5); one biopsy sample of each subject was analyzed without replicates. Finally, only one GPL24557 dataset, GSE140994 (Winkler et al., 2021), which included five normal samples and five NASH samples, was selected as the test set. In addition, six datasets including GSE17470 (Baker et al., 2010), GSE63067 (Frades et al., 2015), GSE35961 (Kita et al., 2012), GSE93819 (Kobori et al., 2017), GSE59492 (Blaya et al., 2016), and GSE33857 (Murakami et al., 2012) were selected as the validation sets (Table 1).

TABLE 1

TABLE 1. Information of selected GEO datasets.

Data normalization and identification of differentially expressed genes (DEGs)

The original files that were downloaded from the GEO databases were preprocessed and normalized by the robust multiarray average (RMA) method by using the GEOquery package in R software (version 4.0.1). The limma R package was used to conduct gene analysis to identify the DEGs of intergroup differences via the “makeContrasts” code. The screening criteria were log2 (fold change) > 2 or < -2 with an adjusted p-value (Q value) < 0.05.

Analysis for heatmap and volcano plot

We performed principal component analysis and drew heatmaps and volcano plots to visualize these DEGs by using the ggplot function and pheatmap package in R software. The online tool BioGPS (http://biogps.org/) was used to assess the tissue-/organ-specific expressed genes. The screening criteria of tissue-specific genes were as follows (1) transcripts that mapped to a specific organ system with an expression value of >10 multiples of the median and (2) second most abundant tissue’s expression were no more than half as high.

Enrichment analysis by GO database

We analyzed the biological process and molecular function according to Gene Ontology (GO) classification to uncover the functions of DEGs and to further predict the biological connections in these genes. The gene list consisting of differentially expressed mRNAs screened was submitted to the DAVID database (https://david.ncifcrf.gov/home.jsp), which is a web-based bioinformatics resource for functional enrichment analysis. The R software cluster profiler package was used for the visualization of the GO enrichment. The top ten GO terms with the most significant p values were selected according to the biological process and molecular function analysis.

Prediction of the PPI network

A protein–protein interaction (PPI) network was predicted based on the DEGs to further explore the functional role of these DEGs by using the online STRING tool (https://strin-g-db.org/) with a filter condition (combined score >0.4) (Szklarczyk et al., 2021). Next, we analyzed the interaction information and drew the PPI network by using the plugin CytoHubba. Minimal Common Oncology Data Elements (MCODE) was used to identify three significant gene clusters and to calculate the cluster scores. Five algorithms, including degree, maximal clique centrality (MCC), maximum neighborhood component (MNC), density of maximum neighborhood component (DMNC), and clustering co-efficient, were utilized to filter the modules of core genes and to obtain the top 15 hub genes. CytoHubba was used to further identify crucial genes in this PPI network as hub genes. The final hub genes were then obtained based on the results of organ-/tissue-specific expression.

Verification of specifically expressed hub genes by four datasets from the GEO database

GSE17470 and GSE63067, which included 16 liver samples from NASH patients and 11 liver samples from control patients, and GSE35961 and GSE93819, which included nine liver samples from NASH mice and nine liver samples from normal mice, were used to verify the specifically expressed hub genes. GraphPad Prism was used to analyze the expression profiles and ROC curves of specifically expressed hub genes.

Prediction of mRNA–miRNA interaction network

We uploaded the specifically expressed hub genes to Network Analyst (http://www.networkanalyst.ca) to explore the potential miRNA–protein interactions. Based on the miRNA target gene prediction, Cytoscape software was used to plot the mRNA–miRNA interaction network and pathway enrichment analysis.

Prediction of ceRNA network

StarBase version 3.0 (http://starbase.sysu.edu.cn/index.php) was used to predict the interaction of lncRNAs and circRNAs with miRNAs (Li et al., 2014). The intersections of the predicted results were used as the target lncRNAs and circRNAs. The ceRNA networks containing lncRNA, miRNA, and mRNA were visualized by the Cytoscape software.

Statistical analysis

Student’s t-test was used to compare the differences between the two groups by IBM SPSS Statistics 25 (SPSS, Inc., Chicago, United States ). All the tests were two-tailed, and p-value less than 0.05 was considered statistically significant. GraphPad Prism 9 was used to perform diagnostic analysis and draw the ROC curves.

Results

Significant damage of LSECs in MCD diet-induced NASH model

Wild-type rats were fed with an MCD diet for 6 weeks to induce NASH. H&E revealed steatosis, ballooning, inflammation, and fibrosis in the livers of MCD-fed rats (Figure 1A). NASH was associated with a remarkable increase in ALT and AST (Figure 1B). Normal LSECs on SEM and TEM imaging were characterized by regularly organized fenestrae and sieve plates on the sinusoidal endothelium (Figures 1C, D). While in the liver of the NASH model, remarkable formation of large defects on sinusoidal endothelium and LSEC lining detachment/deficiency were observed, indicating significant LSEC structural damage (Figures 1C, D).

FIGURE 1

FIGURE 1. LSECs injury in the NASH model. (A) H&E staining was used to evaluate the pathology of the liver in two groups. (B) Alanine aminotransferase (ALT) and aspartate aminotransferase (AST) were tested to evaluate the liver injury. (C,D) Scanning electron microscopy (SEM) and transmission electron microscopy (TEM) showed the ultrastructural features of liver sinusoidal endothelial cells (LSECs) in the normal liver and NASH model. ns, non-significant; ****p < 0.001 by t-test; n = 6 at each time point per group; figures are representative of three experiments.

Identification of DEGs in damaged LSECs via GEO database

As shown in Table 1, a total of seven GEO databases were identified. The dataset GSE140994 contains the gene expression data of LSECs isolated from NASH and NC mouse livers and thus was used as the test set. The datasets, including GSE17470, GSE63067, GSE35961, and GSE93819, contain gene expression profiling of the liver from NASH patients or mice and thus were used as validation sets. GSE59492 and GSE33857 contain non-coding RNA profiling of the liver and serum from NASH patients, respectively, and were used to validate the miRNA expression profile.

Based on the analysis of the dataset GSE140994, a total of 174 DEGs were identified (listed in Supplementary Table S1). As shown in Figure 2A, the principal component analysis indicated a significant difference in the gene expression profile of LSECs between the NASH group and the normal control group. Among the 174 DEGs, there were 156 upregulated genes and 18 downregulated genes (Figure 2B). The heatmap plot was used to visualize these DEGs (Figure 1C).

FIGURE 2

FIGURE 2. Identification of DEGs. (A) Principal component analysis between the normal samples and the NASH samples of LSECs; (B) volcano plot of DEGs between the normal samples and the NASH samples of LSECs. The red plots represent upregulated DEGs, the black plots represent non-significant genes, and the blue plots represent downregulated DEGs. (C) Heatmap of DEGs between the normal samples and the NASH samples of LSECs. Red rectangles represent high expression, and blue rectangles represent low expression.

Identification of the tissue-/organ-specific expressed genes

A total of 84 tissue-/organ-specific expressed genes were identified according to BioGPS (Table 2). We observed that most of these genes were expressed in the immune organ/system (36/84, 42.86%). The digestive system was associated with the second most specifically expressed genes (17/84, 20.24%), followed by the others (10/84, 11.90%), the liver (8/84, 9.52%), the kidney (5/84, 5.95%), and the endocrine system (4/84, 4.76%). Finally, the nervous system and bone were associated with the fewest tissue-/organ-specific expressed genes (2/84, 2.38%) (Table 2).

TABLE 2

TABLE 2. Distribution of tissue- and organ-specific expressed genes identified by BioGPS.

GO enrichment analysis results

The biological process (BP) and molecular function (MF) of these 174 DEGs were analyzed by using the DAVID database. The ten most significantly enriched BP and MF terms were presented in the circle plot (Figure 3A). The outer circle is a scatter plot for each GO term with logFC assigned to genes. The red dots represented upregulated genes, and the blue dots represented downregulated genes. The 10 GO terms included the following: inflammatory response (GO:0006594), immune response (GO:0006955), cell adhesion (GO:0007155), immune system process (GO:0002376), cellular response to interleukin 1 (GO:0071347), chemokine activity (GO:0008009), integrin binding (GO:0005178), cell adhesion molecule binding (GO:0050839), cytokine activity (GO:0005125), and identical protein binding (GO:0042802). A chord plot displayed the distribution of each DEG to the corresponding pathway enrichment (Figure 3B).

FIGURE 3

FIGURE 3. GO (Gene Ontology) terms enrichment analysis. (A) GO circle plot; the inner ring was a bar plot where the bar height showed the significance of the term (−log10 p-value) and the gradual color showed the z-score. The outer ring displayed scatter plots of the expression levels (logFC) for the genes in each term. The distribution of genes in different GO terms was used to predict different annotations. (B) GO chord plot of the relationship between the list of selected genes and their corresponding GO terms, together with the logFC of the genes. The left half of the GO chord displayed the up-expression or down-expression DEGs. The right half represented different GO terms with varied colors. A gene was linked to a certain GO term by the colored bands.

PPI network analysis, MCODE cluster modules, and hub gene identification

An interaction network between the proteins coded by DEGs was predicted by STRING and visualized by Cytoscape (Figure 4A). The PPI network comprised 124 nodes and 526 edges, which indicated each DEG and the interaction between them, respectively. In this network, three cluster modules (Figures 4B–D) were identified by using the MCODE plugin according to the screening criteria. Module 1 (cluster score: 5; with 25 nodes and 60 edges) mainly contains DEGs involved in chemokine activity and expressed in the digestive system. The DEGs in module 2 (cluster score: 5; with 23 nodes and 55 edges) were mainly related to chemokine activity and expressed in the immune system. Module 3 (cluster score: 4.182; with 12 nodes and 23 edges) was mainly related to the signaling molecules of cytokines and chemokines.

FIGURE 4

FIGURE 4. PPI network of DEGs and three clusters extracted by MCODE. (A) Using the Cytoscape software and online database STRING, the protein–protein interaction network coded by DEGs comprised 124 nodes and 526 edges. Each node represents a protein, while each edge represents the association of proteins. Red diamonds represent the upregulated DEGs, and blue hexagons represent the downregulated DEGs. Three cluster modules were extracted by MCODE. Cluster 1 (B) had high cluster score (score: 5, 25 nodes and 60 edges), followed by cluster 2 (C) (score: 5, 23 nodes and 55edges), and cluster 3 (D) (score: 4.182, 12 nodes and 23 edges).

Next, 15 hub genes including CCl4, FN1, CCR2, EPCAM, CDH1, CXCL2, CCL3, KRT19, SPP1, KRT18, C3, CX3CR1, KRT8, ITGAX, and CLDN3 were identified according to the five algorithms of the cytoHubba plugin. The detailed information on these 15 genes was listed in Table 3. As the immune reaction plays an important role in the pathogenesis of NASH and GO enrichment analyses showed that DEGs were mainly enriched in immune-activated and inflammatory pathways, we selected nine hub genes which are specifically expressed in the immune tissues/cells and the digestive system. These nine specifically expressed hub genes could be the most important genes in the PPI network.

TABLE 3

TABLE 3. 15 hub genes identified by five algorithms of cytoHubba.

Prediction of mRNA–miRNA interaction and construction of the co-expressed network

We used the Network Analyst database to predict target miRNAs of hub genes. We obtained 98 target miRNAs for these nine specifically expressed hub genes and 144 mRNA–miRNA connections. A co-expressed network of mRNAs and miRNAs which comprised 107 nodes and 144 edges was plotted by Cytoscape (Figure 5). CXCL2 has the largest number of target miRNAs, while CX3CR1 has only four target miRNAs. Each diamond refers to a specific hub gene, while the spheres connected to the diamond refer to the target miRNAs. The colors of the spheres indicate a varying number of hub genes that are related to the miRNA. Thus, the miRNA hsa-miR-27a-3p with green color is related to the maximum number of hub genes, indicating that it might have the most complex and important biological functions.

FIGURE 5

FIGURE 5. miRNA-mRNA-co-expressed network. A co-expressed network of mRNAs and target miRNAs. The mRNA–miRNA co-expressed network was constructed by Cytoscape including 107 nodes and 144 edges. Every node represents one miRNA/mRNA, and each edge represents the interaction of mRNA and miRNA. mRNAs and miRNAs are indicated with diamond and circle shapes, respectively.

Verification of the nine specifically expressed hub genes

We further validated the expression of the nine specifically expressed hub genes in another four GEO databases as listed in Table 1. When examining the diagnostic performance for the hub genes, we included all the samples from these databases in a diagnostic analysis and calculated the area under the receiver operating characteristic (AUROC). As shown in Figure 6; Supplementary Figure S1, we found an increasing trend of CCL4 expression in the four databases (Figure 6A) with an AUROC value of 0.77 (Figure 6B). Meanwhile, the expression of ITGAX was significantly increased in the four databases (p < 0.05 for all, Figure 6C) with a favorable AUROC value of 0.86 (Figure 6D). On the other hand, for the other seven hub genes, no significantly increased/decreased expression (or trend) was identified in the four databases (Supplementary Figure S1).

FIGURE 6

FIGURE 6. Identification of the CCL4 and ITGAX in four validation databases and investigation of its diagnosis value in NASH. (A) mRNA expression of CCL4 in NASH and normal samples was verified in four databases. (B) ROC curves of the CCL4 in the validating set. (C) mRNA expression of ITGAX in NASH and normal samples were verified in four databases. (D) ROC curves of ITGAX in the validating set. ns, non-significant; *p < 0.05 and ***p < 0.005 by t-test.

Verification of miRNA in the liver and serum of NASH patients

Using the dataset GSE59492, which contains the non-coding RNA profiling of the liver of NASH patients, we finally identified 37 out of 144 miRNAs that were downregulated in the liver of NASH patients (Table 4). Among these 37 miRNAs, there were 23 miRNAs that were also downregulated in the serum of NASH patients (Table 4). Among these miRNA, hsa-miR-378a, hsa-miR-17, and hsa-miR-20a are related to CCL4 in the mRNA–miRNA co-expression network (Figure 5), while hsa-miR-103a and hsa-miR-107 are related to ITGAX (Figure 5).

TABLE 4

TABLE 4. Verification of miRNA in the liver and serum of NASH patients.

Prediction of target non-coding RNAs and ceRNA networks

MiRNAs are known to induce gene silencing and downregulate gene expression by binding to mRNAs. However, its upstream molecules, such as circRNAs and lncRNAs, can affect the function of miRNA by combining miRNA response elements, thus regulating gene expression via the transcript process. This interaction between RNAs is called a ceRNA network. Subsequently, we explored the lncRNAs and circRNAs related to the aforementioned miRNAs by using the Starbase tool (https://starbase.sysu.edu.cn/). Finally, we obtained the networks involving 13 circRNAs, one lncRNA related to the miRNAs of CCL4 (Figure 7A), 15 circRNAs, and one lncRNA related to the miRNAs of ITGAX (Figure 7B). We propose that the lincRNA-circRNA-miRNA-CCL4 and lincRNA-circRNA-miRNA-ITGAX networks might be potential key pathways to regulate the disease of NASH.

FIGURE 7

FIGURE 7. lncRNA/circRNA-miRNA-mRNA-ceRNA network. (A) Network of CCL4-IncRNA-miRNA-mRNA. (B) Network of ITGAX-IncRNA-miRNA-mRNA.

Discussion

In our study, we obtained 174 DEGs of LSECs in NASH. GO enrichment analysis revealed in transcript levels that these DEGs mainly participated in the inflammatory response and cell adhesion activity in NASH. GO pathways that were enriched included inflammatory and immunological pathways involving inflammatory response interaction, cell adhesion, cellular response to the interleukin-1 signaling pathway, and chemokine and cytokine activities. These results suggest that LSECs in NASH have distinct immune activation which is one of the main causes of liver inflammation and fibrosis.

Two immune system-specifically expressed genes, namely, CCL4 and ITGAX, were finally identified. CCL4 and ITGAX were both upregulated in LSECs in NASH with statistical significance (p < 0.05). ROC curves showed that both genes had a certain diagnostic value for both human and animal NASH. CCL4 and ITGAX can be potential biomarkers for early diagnosis of NASH. In addition, we constructed an mRNA–miRNA co-expression network and ceRNA networks to further understand the pathogenesis of NASH at the transcriptome level. The prediction of upstream miRNAs, lncRNAs, and circRNAs led to the identification of a potential mRNA-miRNA-lncRNA/circRNA regulatory network.

ITGAX (integrin alpha x), also known as CD11c, a member of the integrin family, serves as a receptor for the extracellular matrix. ITGAX is a surface marker that can be expressed on dendritic cells, B lymphocytes, T lymphocytes, NK cells, and subsets of monocytes and macrophages (Poltorak and Schraml, 2015). The expression of ITGAX by LSECs has been reported in the previous publication (Knolle et al., 1999). Flow cytometry showed that LSECs expressed ITGAX in vitro, which could be downregulated by treatment of cells with endotoxin (Knolle et al., 1999). This could explain that in situ studies showed the lack of ITGAX expression on LSECs because LSECs are physiologically exposed to endotoxin in portal venous blood. Therefore, under certain circumstances, such as in vitro or hepatopathy, the expression of ITGAX by LSECs might be activated. Recently, accumulating evidence indicates that ITGAX might play an important role in cell death, cytokine processing, and inflammation (Wang et al., 2019). An increased number of ITGAX⁺ cells has been found in the liver of mice under obese conditions (Stefanovic-Racic et al., 2012), which is one of the key triggers for liver inflammation progression from NAFLD to NASH (Younossi et al., 2016). This indicates that ITGAX might be a significant mediator in the pathogenesis of NASH. Our study found that ITGAX was significantly upregulated in the LSECs in NASH. The ROC curve of ITGAX indicated that it had a very strong correlation with NASH (AUROC = 0.86). ITGAX is considered a promising biomarker for the diagnosis of NASH.

CCL4, also known as macrophage inflammatory protein-1β, is a chemo-attractant for varieties of immune cells (Zhao et al., 2021). CCR5, the receptor of CCL4, is expressed by various cells, including T lymphocytes, NK, NKT cells, monocytes, macrophages, and γδ T cells. Quite a few previous reports as well as our previous article have proved that LSEC is capable of expressing inflammatory mediators including cytokines and chemokines (Wang et al., 2021). CCL4 has also been found to be expressed in vascular endothelial cells (Barker et al., 2008). It is found that inflammatory factors, such as IL-1beta, can mediate the upregulation of CCL4 production in the hepatic cells, thus contributing to the continuous recruitment of activating immunocytes to the liver (Zhang et al., 2003). The increased CCL4 is related to critical biological processes of inflammation in the liver. In addition, miRNAs could upregulate CCL4 expression in immune cells (Cheng et al., 2015). In our study, we also found several circRNAs and lncRNAs that were related to the regulation of CCL4 expression in LSECs. Therefore, further investigation might as well focus on illustrating the mechanism of the CCL4-mRNA-miRNA-lncRNA/circRNA regulatory network in NASH. LncRNAs have emerged as significant factors in almost all processes of gene function and regulation. In addition, it has been reported that the Mir503HG is involved in tumor metastasis in hepatocellular carcinoma (Wang et al., 2018). Therefore, MIR503HG-hsa-miR-107/hsa-miR-103a/-ITGAX might be a potential ceRNA regulatory pathway to regulate the disease progression of NASH.

There are limitations to this study. The major disadvantage of the MCD diet model is the absence of insulin resistance and metabolic syndrome. In an MCD diet model, loss of body weight is the converse of obesity-related NASH. Meanwhile, in the MCD diet model, initial fatty changes can be observed in the mice liver after 2 weeks of MCD diet followed by progressive deterioration to evident NASH 6 weeks after the beginning of treatment. Thus, initial simple steatosis can progress into NASH within 6 weeks in the MCD diet model. The liver pathology examination also revealed that the MCD diet model recapitulates the major characteristics of human NASH, including steatosis, ballooning degeneration, inflammation, and fibrosis. Therefore, the MCD diet model performs well in mimicking the changes in the intrahepatic environment, especially the liver-specific cells such as LSECs. It would be better that the changes of LSECs can be further validated in other models such as high-fat diet-induced NAFLD/NASH.

Conclusion

In this study, we systematically analyzed the differentially expressed genes of LSECs in NASH. We identified two immune system-specifically expressed genes, CCL4 and ITAGX, expressed by LSECs, as potential key players for the pathogenesis of disease development in NASH at the transcriptome level. Our findings put forward promising therapeutic and prognostic values of the genes expressed by LSECs as well as the ceRNA regulatory networks in NASH.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

Ethics statement

The animal study was reviewed and approved by The Ethics Committee of Peking University People’s Hospital.

Author contributions

(I) Conception and design: YulL and YW. (II) Administrative support: YulL and YW. (III) Provision of study materials or patients: YZ, YimL, and YunL. (IV) Collection and assembly of data: YZ, YimL, YunL, and YW. (V) Data analysis and interpretation: YZ, YimL, YulL, and YW. (VI) Manuscript writing: All authors. (VII) Final approval of manuscript: All authors.

Funding

This work is supported by the National Natural Science Foundation of China (Grant numbers 82070539, 81873549, and 82100579). The sponsors had no involvement in the collection, analysis, and interpretation of data; in the writing of the report; and in the decision to submit the article for publication.

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/fcell.2022.946566/full#supplementary-material

Abbreviations

NASH - nonalcoholic steatohepatitis; NAFLD - nonalcoholic fatty liver disease; LSECs - liver sinusoidal endothelial cells; NK - natural killer; LPS - lipopolysaccharide; ceRNA - competitive endogenous RNA; MCD - methionine–choline-deficient diet; ALT - alanine aminotransferase; AST - aspartate aminotransferase; SEM - scanning electron microscopy; TEM - transmission electron microscopy; GO - Gene Ontology; DEG - differentially expressed genes; PPI - protein–protein interaction; MCODE - Minimal Common Oncology Data Elements; BP - biological process; MF - molecular function; ITGAX - integrin alpha x.

References

Baker, S. S., Baker, R. D., Liu, W., Nowak, N. J., and Zhu, L. (2010). Role of alcohol metabolism in non-alcoholic steatohepatitis. PloS one 5 (3), e9570. PubMed PMID: 20221393; PubMed Central PMCID: PMC2833196. doi:10.1371/journal.pone.0009570

PubMed Abstract | CrossRef Full Text | Google Scholar

Barker, K. S., Park, H., Phan, Q. T., Xu, L., Homayouni, R., Rogers, P. D., et al. (2008). Transcriptome profile of the vascular endothelial cell response to Candida albicans. J. Infect. Dis. 198 (2), 193–202. PubMed PMID: 18500935. doi:10.1086/589516

PubMed Abstract | CrossRef Full Text | Google Scholar

Blaya, D., Coll, M., Rodrigo-Torres, D., Vila-Casadesus, M., Altamirano, J., Llopis, M., et al. (2016). Integrative microRNA profiling in alcoholic hepatitis reveals a role for microRNA-182 in liver injury and inflammation. Gut 65 (9), 1535–1545. PubMed PMID: 27196584. doi:10.1136/gutjnl-2015-311314

PubMed Abstract | CrossRef Full Text | Google Scholar

Cai, J., Zhang, X. J., and Li, H. (2019). The role of innate immune cells in nonalcoholic steatohepatitis. Hepatology 70 (3), 1026–1037. PubMed PMID: 30653691. doi:10.1002/hep.30506

PubMed Abstract | CrossRef Full Text | Google Scholar

Carambia, A., Freund, B., Schwinge, D., Heine, M., Laschtowitz, A., Huber, S., et al. (2014). TGF-β-dependent induction of CD4⁺CD25⁺Foxp3⁺ Tregs by liver sinusoidal endothelial cells. J. Hepatol. 61 (3), 594–599. PubMed PMID: 24798620. doi:10.1016/j.jhep.2014.04.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, N. L., Chen, X., Kim, J., Shi, A. H., Nguyen, C., Wersto, R., et al. (2015). MicroRNA-125b modulates inflammatory chemokine CCL4 expression in immune cells and its reduction causes CCL4 increase with age. Aging Cell 14 (2), 200–208. PubMed PMID: 25620312; PubMed Central PMCID: PMC4364832. doi:10.1111/acel.12294

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, Q., Chen, X., Wu, H., and Du, Y. (2021). Three hematologic/immune system-specific expressed genes are considered as the potential biomarkers for the diagnosis of early rheumatoid arthritis through bioinformatics analysis. J. Transl. Med. 19 (1), 18. PubMed PMID: 33407587; PubMed Central PMCID: PMC7789535. doi:10.1186/s12967-020-02689-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Demircioglu, D., Cukuroglu, E., Kindermans, M., Nandi, T., Calabrese, C., Fonseca, N. A., et al. (2019). A pan-cancer transcriptome analysis reveals pervasive regulation through alternative promoters. Cell 178 (6), 1465–1477. e17PubMed PMID: 31491388. doi:10.1016/j.cell.2019.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Dufour, J. F., Caussy, C., and Loomba, R. (2020). Combination therapy for non-alcoholic steatohepatitis: Rationale, opportunities and challenges. Gut 69 (10), 1877–1884. PubMed PMID: 32381514; PubMed Central PMCID: PMC7497577. doi:10.1136/gutjnl-2019-319104

PubMed Abstract | CrossRef Full Text | Google Scholar

Frades, I., Andreasson, E., Mato, J. M., Alexandersson, E., Matthiesen, R., and Martinez-Chantar, M. L. (2015). Integrative genomic signatures of hepatocellular carcinoma derived from nonalcoholic Fatty liver disease. PloS one 10 (5), e0124544. PubMed PMID: 25993042; PubMed Central PMCID: PMC4439034. doi:10.1371/journal.pone.0124544

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujimoto, D., Kuwabara, T., Hata, Y., Umemoto, S., Kanki, T., Nishiguchi, Y., et al. (2020). Suppressed ER-associated degradation by intraglomerular cross talk between mesangial cells and podocytes causes podocyte injury in diabetic kidney disease. FASEB J. official Publ. Fed. Am. Soc. Exp. Biol. 34 (11), 15577–15590. PubMed PMID: 32996639. doi:10.1096/fj.202000078RR

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammoutene, A., and Rautou, P. E. (2019). Role of liver sinusoidal endothelial cells in non-alcoholic fatty liver disease. J. Hepatol. 70 (6), 1278–1291. PubMed PMID: 30797053. doi:10.1016/j.jhep.2019.02.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Hilscher, M. B., Sehrawat, T., Arab, J. P., Zeng, Z., Gao, J., Liu, M., et al. (2019). Mechanical stretch increases expression of CXCL1 in liver sinusoidal endothelial cells to recruit neutrophils, generate sinusoidal microthombi, and promote portal hypertension. Gastroenterology 157 (1), 193–209. e9PubMed PMID: 30872106; PubMed Central PMCID: PMC6581607. doi:10.1053/j.gastro.2019.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Kita, Y., Takamura, T., Misu, H., Ota, T., Kurita, S., Takeshita, Y., et al. (2012). Metformin prevents and reverses inflammation in a non-diabetic mouse model of nonalcoholic steatohepatitis. PloS one 7 (9), e43056. PubMed PMID: 23028442; PubMed Central PMCID: PMC3445596. doi:10.1371/journal.pone.0043056

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleiner, D. E., Brunt, E. M., Wilson, L. A., Behling, C., Guy, C., Contos, M., et al. (2019). Association of histologic disease activity with progression of nonalcoholic fatty liver disease. JAMA Netw. Open 2 (10), e1912565. PubMed PMID: 31584681; PubMed Central PMCID: PMC6784786. doi:10.1001/jamanetworkopen.2019.12565

PubMed Abstract | CrossRef Full Text | Google Scholar

Knolle, P. A., Schmitt, E., Jin, S., Germann, T., Duchmann, R., Hegenbarth, S., et al. (1999). Induction of cytokine production in naive CD4(+) T cells by antigen-presenting murine liver sinusoidal endothelial cells but failure to induce differentiation toward Th1 cells. Gastroenterology 116 (6), 1428–1440. PubMed PMID: 10348827. doi:10.1016/s0016-5085(99)70508-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobori, M., Takahashi, Y., Sakurai, M., Ni, Y., Chen, G., Nagashimada, M., et al. (2017). Hepatic transcriptome profiles of mice with diet-induced nonalcoholic steatohepatitis treated with astaxanthin and vitamin E. Int. J. Mol. Sci. 18 (3), E593. PubMed PMID: 28282876; PubMed Central PMCID: PMC5372609. doi:10.3390/ijms18030593

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishnan, M. S., Rajan, Kd A., Park, J., Arjunan, V., Garcia Marques, F. J., Bermudez, A., et al. (2021). Genomic analysis of vascular invasion in HCC reveals molecular drivers and predictive biomarkers. Hepatology 73 (6), 2342–2360. PubMed PMID: 33140851; PubMed Central PMCID: PMC8115767. doi:10.1002/hep.31614

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J. H., Liu, S., Zhou, H., Qu, L. H., and Yang, J. H. (2014). starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 42, D92–D97. PubMed PMID: 24297251; PubMed Central PMCID: PMC3964941. doi:10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Murakami, Y., Toyoda, H., Tanahashi, T., Tanaka, J., Kumada, T., Yoshioka, Y., et al. (2012). Comprehensive miRNA expression analysis in peripheral blood can diagnose liver disease. PloS one 7 (10), e48366. PubMed PMID: 23152743; PubMed Central PMCID: PMC3485241. doi:10.1371/journal.pone.0048366

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng-Nguyen, D., Van Nguyen, T., Van Nguyen, T., Nguyen, H. Q., and Nguyen, V. T. (2021). Prevalence and risk factors of Taenia hydatigena in dogs, pigs, and cattle in the Central Highlands of Vietnam. Parasitol. Res. 120 (9), 3245–3253. PubMed PMID: 34386855. doi:10.1007/s00436-021-07276-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Poltorak, M. P., and Schraml, B. U. (2015). Fate mapping of dendritic cells. Front. Immunol. 6, 199. PubMed PMID: 25999945; PubMed Central PMCID: PMC4418393. doi:10.3389/fimmu.2015.00199

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefanovic-Racic, M., Yang, X., Turner, M. S., Mantell, B. S., Stolz, D. B., Sumpter, T. L., et al. (2012). Dendritic cells promote macrophage infiltration and comprise a substantial proportion of obesity-associated increases in CD11c+ cells in adipose tissue and liver. Diabetes 61 (9), 2330–2339. PubMed PMID: 22851575; PubMed Central PMCID: PMC3425427. doi:10.2337/db11-1523

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Gable, A. L., Nastou, K. C., Lyon, D., Kirsch, R., Pyysalo, S., et al. (2021). The STRING database in 2021: Customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 49 (D1), D605–D612. PubMed PMID: 33237311; PubMed Central PMCID: PMC7779004. doi:10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Liang, L., Dong, Q., Huan, L., He, J., Li, B., et al. (2018). Long noncoding RNA miR503HG, a prognostic indicator, inhibits tumor metastasis by regulating the HNRNPA2B1/NF-κB pathway in hepatocellular carcinoma. Theranostics 8 (10), 2814–2829. PubMed PMID: 29774077; PubMed Central PMCID: PMC5957011. doi:10.7150/thno.23012

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Yang, L., Liang, F., Chen, Y., and Yang, G. (2019). Integrin alpha x stimulates cancer angiogenesis through PI3K/Akt signaling-mediated VEGFR2/VEGF-A overexpression in blood vessel endothelial cells. J. Cell. Biochem. 120 (2), 1807–1818. PubMed PMID: 30873824. doi:10.1002/jcb.27480

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., and Liu, Y. (2021). Gut-liver-axis: Barrier function of liver sinusoidal endothelial cell. J. Gastroenterol. Hepatol. 36 (10), 2706–2714. PubMed PMID: 33811372. doi:10.1111/jgh.15512

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., and Liu, Y. (2021). Neutrophil-induced liver injury and interactions between neutrophils and liver sinusoidal endothelial cells. Inflammation 44 (4), 1246–1262. PubMed PMID: 33649876. doi:10.1007/s10753-021-01442-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Zhang, Y., Liu, Y., Xu, J., and Liu, Y. (2021). Gut-liver Axis: Liver sinusoidal endothelial cells function as the hepatic barrier in colitis-induced liver injury. Front. Cell Dev. Biol. 9, 702890. PubMed PMID: 34336855; PubMed Central PMCID: PMC8322652. doi:10.3389/fcell.2021.702890

PubMed Abstract | CrossRef Full Text | Google Scholar

Winkler, M., Staniczek, T., Kurschner, S. W., Schmid, C. D., Schonhaber, H., Cordero, J., et al. (2021). Endothelial GATA4 controls liver fibrosis and regeneration by preventing a pathogenic switch in angiocrine signaling. J. Hepatol. 74 (2), 380–393. PubMed PMID: 32916216. doi:10.1016/j.jhep.2020.08.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Wisse, E., De Wilde, A., and De Zanger, R. (1984). “Perfusion fixation of human and rat liver tissue for light and electron microscopy: A review and assessment of existing methods with special emphasis on sinusoidal cells and microcirculation,” in The science of biological specimen preparation for microscopy and microanalysis (Switzerland: Springer), 31–38.

Google Scholar

Younossi, Z. M., Koenig, A. B., Abdelatif, D., Fazel, Y., Henry, L., and Wymer, M. (2016). Global epidemiology of nonalcoholic fatty liver disease-Meta-analytic assessment of prevalence, incidence, and outcomes. Hepatology 64 (1), 73–84. PubMed PMID: 26707365. doi:10.1002/hep.28431

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, T., Guo, C. J., Li, Y., Douglas, S. D., Qi, X. X., Song, L., et al. (2003). Interleukin-1beta induces macrophage inflammatory protein-1beta expression in human hepatocytes. Cell. Immunol. 226 (1), 45–53. PubMed PMID: 14746807; PubMed Central PMCID: PMC4016814. doi:10.1016/j.cellimm.2003.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, N., Dang, H., Ma, L., Martin, S. P., Forgues, M., Ylaya, K., et al. (2021). Intratumoral γδ T-cell infiltrates, chemokine (C-C motif) ligand 4/chemokine (C-C motif) ligand 5 protein expression and survival in patients with hepatocellular carcinoma. Hepatology 73 (3), 1045–1060. PubMed PMID: 32502310. doi:10.1002/hep.31412

PubMed Abstract | CrossRef Full Text | Google Scholar