Identification of miRNA-target gene regulatory networks in liver fibrosis based on bioinformatics analysis

Background Liver cirrhosis is one of the leading causes of death worldwide. MicroRNAs (miRNAs) can regulate liver fibrosis, but the underlying mechanisms are not fully understood, and the interactions between miRNAs and mRNAs are not clearly elucidated. Methods miRNA and mRNA expression arrays of cirrhotic samples and control samples were acquired from the Gene Expression Omnibus database. miRNA-mRNA integrated analysis, functional enrichment analysis and protein-protein interaction (PPI) network construction were performed to identify differentially expressed miRNAs (DEMs) and mRNAs (DEGs), miRNA-mRNA interaction networks, enriched pathways and hub genes. Finally, the results were validated with in vitro cell models. Results By bioinformatics analysis, we identified 13 DEMs between cirrhotic samples and control samples. Among these DEMs, six upregulated (hsa-miR-146b-5p, hsa-miR-150-5p, hsa-miR-224-3p, hsa-miR-3135b, hsa-miR-3195, and hsa-miR-4725-3p) and seven downregulated (hsa-miR-1234-3p, hsa-miR-30b-3p, hsa-miR-3162-3p, hsa-miR-548aj-3p, hsa-miR-548x-3p, hsa-miR-548z, and hsa-miR-890) miRNAs were further validated in activated LX2 cells. miRNA-mRNA interaction networks revealed a total of 361 miRNA-mRNA pairs between 13 miRNAs and 245 corresponding target genes. Moreover, PPI network analysis revealed the top 20 hub genes, including COL1A1, FBN1 and TIMP3, which were involved in extracellular matrix (ECM) organization; CCL5, CXCL9, CXCL12, LCK and CD24, which participated in the immune response; and CDH1, PECAM1, SELL and CAV1, which regulated cell adhesion. Functional enrichment analysis of all DEGs as well as hub genes showed similar results, as ECM-associated pathways, cell surface interaction and adhesion, and immune response were significantly enriched in both analyses. Conclusions We identified 13 differentially expressed miRNAs as potential biomarkers of liver cirrhosis. Moreover, we identified 361 regulatory pairs of miRNA-mRNA and 20 hub genes in liver cirrhosis, most of which were involved in collagen and ECM components, immune response, and cell adhesion. These results would provide novel mechanistic insights into the pathogenesis of liver cirrhosis and identify candidate targets for its treatment.


INTRODUCTION
Liver cirrhosis is a serious public health problem in the world, which causes approximately 1 million human deaths each year (Asrani et al., 2019). The most common primary etiologies for cirrhosis include hepatitis B virus (HBV), hepatitis C virus (HCV), as well as nonalcoholic fatty liver disease (NAFLD) and alcoholic liver disease (ALD) (Asrani et al., 2019). In the context of liver injury and inflammation, excessive amounts of extracellular matrix (ECM) are deposited in the liver and play a central role in the pathogenesis of fibrosis and cirrhosis (Friedman, 2003). This process is regulated by various cells and molecular pathways, such as the activation of hepatic stellate cells (HSCs) and portal fibroblasts (Iwaisako et al., 2014;Nishio et al., 2019), infiltration of immune cells (Gehrke et al., 2018;Krenkel & Tacke, 2017), and dedifferentiation of liver sinusoidal epithelial cells (Ding et al., 2014). However, the mechanisms underlying the development of liver fibrosis are not fully understood. Currently, no treatment can effectively reverse decompensated liver cirrhosis. Therefore, studies further elucidating the pathophysiological mechanisms of the disease are warranted to provide hints for treatment of human liver cirrhosis.
MicroRNAs (miRNAs) are a class of small, noncoding RNAs that can bind to messenger RNAs (mRNAs) and repress gene expression (Ambros, 2004). miRNAs play pivotal roles in regulating diverse biological and pathological processes, including cell proliferation (Shenoy & Blelloch, 2014), metabolism (Agbu & Carthew, 2021), inflammation (Naqvi et al., in press) and disease progression, especially chronic liver diseases (Szabo & Bala, 2013;Wang et al., 2021). miRNAs specifically enriched in the liver can regulate hepatic glucose and lipid metabolism, and their expression is correlated with the disease severity of nonalcoholic steatohepatitis (NASH) (Miyaaki et al., 2014;Yamada et al., 2013). Regarding the effect of miRNAs on liver fibrosis, previous studies have yielded inconsistent results, as some showed that miRNAs were profibrogenic, while others showed that miRNAs exerted antifibrogenic effects. For example, Zhang et al. found that inhibition of miR-21, a miRNA enriched in fibrotic tissues, dramatically reduced HSC activation and ameliorated liver fibrosis by inhibiting the transforming growth factor-β (TGF-β) pathway (Zhang et al., 2013). Moreover, treatment with miR-29a, the miRNA most highly expressed in HSCs, significantly attenuated liver fibrosis in animal models by downregulating ECM-related gene expression (Matsumoto et al., 2016;Roderburg et al., 2011). In contrast, Caviglia et al. found that neither the specific deletion of miR-21 in HSCs nor the pan-deletion of miRNAs in HSCs exerted significant effects on HSC activation and liver fibrosis (Caviglia et al., 2018). As shown, the regulatory effect of miRNAs on liver fibrosis is still under debate. Moreover, little is known about the downstream targets and biological pathways involved in the regulation of liver fibrosis by miRNAs.
In this study, we analyzed two microarray datasets to identify differentially expressed miRNAs and mRNAs between cirrhotic patients and healthy controls. Moreover, miRNA-mRNA integrated analysis, functional enrichment analysis and protein-protein interaction (PPI) network construction were performed to better understand the downstream effects of miRNA-mRNA regulation. Finally, the results were validated with in vitro cell models.
The mRNA microarray data of GSE14323 were also extracted from the GEO database and were generated using an Affymetrix Human Genome U133A Array (GPL96) and Affymetrix Human Genome U133A 2.0 Array (GPL571). A total of 36 liver samples were included in this dataset, 17 of which were from cirrhotic patients (Cirrhosis) whereas 19 were from normal controls (Normal).

Identification of differentially expressed miRNAs and mRNAs
The raw data were preprocessed by background correction, gene symbol transformation, and normalization through R programming. Then, the Limma R package (v3.38.3) was used to identify the differentially expressed miRNAs (DEMs) and mRNAs (DEGs) between cirrhotic patients and controls with an unpaired t -test. With the thresholds of P < 0.05 and |fold change| ≥ 1.5, DEMs were identified between ALD-CH and Normal, HCV-CH and Normal, and NASH-CH and Normal. A Venn diagram was generated using the online tool jvenn (http://www.bioinformatics.com.cn/static/others/jvenn/) to overlap the DEMs in ALD-CH, HCV-CH, and NASH-CH. In addition, mRNAs with P < 0.05 and |fold change|≥ 2 were selected as the significant DEGs between cirrhotic samples and normal samples. Additionally, volcano plots and heatmaps were generated to visualize the significant gene expression changes using GraphPad Prism 8.4.0 and TB (Toolbox for Biologists) tools software (v1.082) (Chen et al., 2020a), respectively.

miRNA-mRNA integrated analysis
The target genes of DEMs were predicted using the online tool miRWalk 3.0 (http: //mirwalk.umm.uni-heidelberg.de/) (Sticht et al., 2018). The overlapping genes between the predicted DEM targets and DEGs were selected for integrated analysis. Since miRNAs generally suppress the expression of target mRNAs, only those differentially expressed target genes that were negatively correlated with the specific miRNA were subjected to subsequent analysis ( Fig. 1). Then, the miRNA-mRNA interaction networks in liver cirrhosis were constructed and visualized using Cytoscape software (v3.8.2).

Functional enrichment analysis
Using Metascape (https://metascape.org/) (Zhou et al., 2019), GO analysis was performed to annotate the biological process (BP), cellular component (CC) and molecular function (MF) of the dysregulated genes, and KEGG and Reactome pathway enrichment were used to identify the critical signal pathways of the DEGs. Enrichment analysis was plotted by http://www.bioinformatics.com.cn and https://cloud.oebiotech.cn, two online platforms for data analysis and visualization. Any GO terms and pathways with a P value less than 0.05 were considered significantly enriched.

PPI network construction and hub genes identification
The differentially expressed target genes were mapped to the STRING database (https://string-db.org/) (Szklarczyk et al., 2019) to screen PPIs, and the results were visualized by Cytoscape software as a network structure. Each node represented a protein encoded by a DEG, and the edges between nodes represented interactions between proteins. The cytoHubba plug-in of Cytoscape was used to calculate the degree rank of the hub genes, and the top 20 were identified for further enrichment analysis and PPI subnetwork construction.

Cell culture and treatment
The human HSC line LX2 was obtained from Procell Life Science and Technology Co., Ltd.
(Wuhan, China). The cells were cultured in Dulbecco's modified Eagle's medium (DMEM, HyClone, USA) supplemented with 10% fetal bovine serum (FBS, Biological Industries, USA), and 100 U/mL penicillin and 100 µg/mL streptomycin (HyClone) in a humidified atmosphere at 37 • C with 5% CO 2 in air. After 16 h of serum starvation, LX2 cells were treated with 5 ng/mL TGF-β1 (PeproTech, USA) for 24 h and then harvested for gene and protein analysis.

Western blotting
Whole proteins from cells were extracted using a protein extraction kit (Nanjing KeyGen Biotech Co., Ltd, Nanjing, China). SDS-PAGE (10%) was used to separate equal amounts of protein from each sample, and then samples were transferred to PVDF membranes (Merck Millipore, USA). The membranes were blocked with 5% nonfat milk and incubated with anti α-SMA (1:1000, Cat# NB600-536, Novus Biologicals) and anti-GAPDH (1:2000, Cat# sc-166574, Santa Cruz Biotechnology, USA) antibodies overnight at 4 • C. After the PVDF membranes were washed and incubated with HRP-conjugated secondary antibodies (1:20000, Cat# ZB-2305, ZSGB-BIO) at room temperature for 2 h, the protein bands were visualized by chemiluminescence using Western Blotting Luminol Reagent (Santa Cruz Biotechnology), and densitometric analyses were performed using Quantity One software (v4.6.2). The protein levels were normalized against GAPDH and were shown as fold changes relative to the control group.

Quantitative real-time PCR (qRT-PCR)
Total RNA was extracted from the cell lysates using a miRNA Isolation Kit (Foregene, Chengdu, China) according to the manufacturer's instructions. Then, equal amounts of 1 µg RNA were reverse transcribed into cDNA using the miRcute Plus miRNA First-Strand cDNA Synthesis Kit (Tiangen, Beijing, China). qRT-PCR was performed in triplicate using SYBR Green qPCR Master Mix (Bimake, USA) on the CFX96 Real-Time PCR Detection System (Bio-Rad). The relative miRNA expression levels were normalized to the U6 snRNA expression level using the 2 − CT method, and the results were shown as fold changes relative to the control group. The primer sequences were listed in Table 1.

Statistical analysis
All data were expressed as mean ± standard deviation and were analyzed using SPSS 19.0 software (SPSS, USA). Student's t -test was performed for comparisons between two groups, and a P value < 0.05 was considered statistically significant.
On the other hand, based on the screening criteria of P < 0.05 and |fold change| ≥ 2, 672 mRNAs (482 upregulated and 190 downregulated) were shown to be differentially expressed between cirrhotic samples and control samples in the GSE14323 dataset (Fig. 3A). The volcano plots and heatmap of the DEGs were shown in Figs. 3A and 3B, respectively.

miRNA-mRNA interaction networks
Following gene expression analysis, 10,604 potential target genes of DEMs were predicted using miRWalk. By considering the overlapping genes, we identified 362 overlapping genes between the predicted target genes and DEGs, which were used for further correlation analysis (Fig. 3C). After filtering out those pairs in which miRNA and mRNA expression shared the same change trend, a total of 361 miRNA-mRNA pairs between 13 miRNAs and 245 corresponding target genes were identified. In these negatively correlating interactions visualized with Cytoscape, 103 regulatory pairs were composed of 6 upregulated miRNAs and 73 downregulated mRNAs (Fig. 4A, Table 3), and 258 pairs were composed of 7 downregulated miRNAs and 172 upregulated mRNAs (Fig. 4B, Table 3). Notably, a single miRNA might interact with different genes, while one gene might also become the   shared target of different miRNAs. Except for hsa-miR-548x-3p, the other DEMs were all predicted to interact with multiple target genes. Additionally, CYTIP, ENAH, FAM129A, FOSB, IFI44L, MITF, PFKFB3, PTGIS, SLA and TMEM100 were identified as targets of hsa-miR-1234-3p, hsa-miR-30b-3p and hsa-miR-3162-3p. These findings revealed a complex regulatory network between miRNAs and mRNAs in the process of liver cirrhosis.

GO enrichment and pathway analyses
To understand the biological function and key pathways of the DEGs, GO annotation, KEGG pathway and Reactome pathway analyses of 245 dysregulated genes from negative regulatory interactions were performed utilizing the Metascape online tool. The top 10 significantly enriched GO terms, including BP, CC and MF, were presented in Figs. 5A-5C. The significantly enriched entries for BP were regulation of cell adhesion, ECM organization and response to mechanical stimulus (Fig. 5A). Moreover, the ECM, side of membrane and endoplasmic reticulum lumen accounted for the majority of the CC terms (Fig. 5B). The most enriched MFs were functions in heparin binding, glycosaminoglycan binding and ECM structural constituent (Fig. 5C). Moreover, the DEGs were also mapped to the KEGG and Reactome pathway databases. As shown in Figs. 5D-5E, immune response, protein and lipid metabolism, ECM organization, cell surface interaction, and signal transduction were among the top 20 most statistically significant enriched pathways.  SPON1, SPRY1, SWAP70, TAGLN, TCF4, THBD, THBS1,  THBS2, TIMP3, TM4SF4, TMEM100, TNFAIP3, TPM1,  TRAF5, TUBB6, TXNIP TBC1D9, TCF4, THBD,  TIMP3, TIPARP, TMEM100, TNFAIP3, TNFRSF21, TPM1,  TUSC3, UGCG, VCAN, WIPF1 (continued on next page)

PPI network and hub genes analyses
To further elucidate the interactions among 245 DEGs, we performed PPI network analysis using logical data originating from the STRING database. After removing the unconnected nodes, a PPI network consisting of 180 nodes and 482 edges was constructed and visualized by Cytoscape (Fig. 6). The top 20 hub genes, which were defined as genes that played essential roles in the network, were distinguished according to the degree calculated by cytoHubba of Cytoscape (Table 4). The PPI subnetwork analysis revealed an average combined score of 0.64, indicating medium-high-degree interactions between these hub genes (low confidence: scores < 0.4; medium: 0.4 to 0.7; high: > 0.7, Fig. 7A) (Von Mering et al., 2005). Among the 20 hub genes, COL1A1, FBN1, and TIMP3 were involved in ECM organization, CCL5, CXCL9, CXCL12, LCK and CD24 participated in the immune response, and CDH1, PECAM1, SELL and CAV1 regulated cell adhesion. In support of this, GO annotation analysis showed that hub genes were closely associated with leukocyte migration, ECM, and cell adhesion molecule binding (Fig. 7B). Additionally, pathway enrichment analysis was also performed based on the Reactome pathway database. In total, 15 Reactome pathways, involving ECM organization, chemokines, cell surface interactions and interleukin signaling, were significantly enriched (Fig. 7C).

DISCUSSION
In this study, we identified 13 differentially expressed miRNAs by bioinformatics analysis. Seven of these miRNAs were further confirmed by in vitro experiments. In addition, 361 pairs of miRNA-mRNA interactions and 20 hub genes among 245 differentially expressed genes were identified, with multiple pathways involved. Among these miRNAs, only some have been studied in organ fibrosis. miR-146b-5p has been mainly focused on the heart, and it was found to activate fibroblast proliferation, migration and fibroblast-to-myofibroblast transition, and its silencing could improve cardiac remodeling (Liao et al., 2021). miR-3135b was found to be elevated in hypertrophic scars (Zhang et al., 2020). In addition, miR-150-5p expression was generally positively correlated with fibrosis severity in rats with carbon tetrachloride-induced liver fibrosis (Chen et al., 2020b). On the other hand, in chronic HBV-infected patients, miR-3162-3p expression was higher in the serum of individuals with significant liver fibrosis than those with nonsignificant liver fibrosis (Zhang et al., 2021). In addition to these aforementioned miRNAs and somewhat elusive information, 11 out of 13 miRNAs were shown to be involved in liver cirrhosis for the first time with bioinformatics analysis, and the expression of 7 miRNAs was proved to undergo significant changes in LX2 cells challenged with TGF-β1 in vitro. Moreover, it was also concluded in this study that miR-150-5p expression was elevated in both cirrhotic tissues and activated LX2 cells, adding more evidence to the debate about whether miR-150-5p expression increases during liver fibrosis (Cabantous et al., 2017;Cai et al., 2018;Chen et al., 2020b). Liver cirrhosis is the end stage of fibrosis, and ECM and collagen are excessively produced and accumulated in the liver throughout fibrosis (Kisseleva & Brenner, 2021). Our GO enrichment, pathway and hub gene analyses undoubtedly showed that these composition shifts were indispensable in liver cirrhosis. Liver injury is the initiating step of liver cirrhosis, and resident macrophages can detect the injury and then release proinflammatory factors to attract other immune cells (Guillot & Tacke, 2019). In addition, components of the ECM, such as hyaluronan, are implicated in immune regulation. The secretion of hyaluronan has been closely associated with maintaining the myofibroblast phenotype (Meran et al., 2007), and high molecular weight hyaluronan promotes immune tolerance by enhancing the numbers of CD4 + CD25 + regulatory T cells (Bollyky et al., 2009). The most significantly enriched pathway in Reactome pathway analysis, signaling by receptor tyrosine kinases, was also recently reported to be associated with liver fibrosis. Fibroblast growth factor receptor, a major family of receptor tyrosine kinases, plays a key role in the progression and resolution of liver fibrosis by interacting with HSCs, hepatocytes and immune cells (Seitz & Hellerbrand, 2021). In addition, knockdown of Gab1, an important protein in the receptor tyrosine kinase pathway, significantly exacerbated liver fibrosis in animal models (Mizutani et al., 2021). According to KEGG pathway analysis, phagosome was most highly enriched pathway, which was consistent with other studies that revealed that autophagy and LC3-associated phagocytosis were closely related to the pathogenesis of liver fibrosis (Li et al., 2020;Wan et al., 2020). In this study, 20 hub genes in the PPI network were identified. Among them, CCL5, CXCL9, CXCL12, LCK and CD24, which are thought to participate in the immune response, were identified as the most affected genes, suggesting a disruption of the immune response in liver cirrhosis. Local inflammation in the liver could cause immune cells to migrate to the foci with the assistance of cell adhesion molecules. Moreover, activated myofibroblastic The data were presented as mean ± standard deviation. n = 3/group, * P < 0.05. α-SMA, α-smooth muscle actin; DAPI, 4 , 6-diamino-2-phenylindole; DEMs, differentially expressed miRNAs; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; HSCs, hepatic stellate cells; qRT-PCR, quantitative real-time PCR; TGF-β1, transforming growth factor-β1.
Full-size DOI: 10.7717/peerj.11910/ fig-8 HSCs gradually transform their surrounding ECM into a high-density matrix with increased stiffness; therefore, cell-matrix interactions would be adapted to these changes (Hintermann & Christen, 2019). Our results showed that the expression of cell adhesion molecules, such as CDH1, PECAM1, SELL and CAV1, which are canonical molecules in cell adhesion, was evidently changed in liver cirrhosis. Of note, platelet-associated pathways were significantly enriched in Reactome pathway analysis, coinciding with former studies that showed that platelet-derived TGF-β is closely involved in the pathogenesis of liver fibrosis (Ahamed & Laurence, 2017;Karolczak & Watala, 2021). Liver sinusoidal epithelial cells can help internalize platelets, thus triggering the profibrogenic effect of platelets (Poisson et al., 2017), which was also reflected in Reactome analysis, as cell surface interactions at the vascular wall was enriched.
There are still several limitations in the present study. The number of samples from GSE59492 and GSE14323 was not large enough, inevitably including some bias during analysis. The transcriptome data being analyzed were from liver biopsies, the differential expression of these miRNAs in cirrhosis versus healthy liver might be accounted for by differential expression in hepatocytes, HSCs, endothelial cells, and Kupffer cells, etc., or through more complex interplays within liver tissues, which could be part of the reason why we did not detect differential expression of all identified miRNAs in LX2 cells. Although HSCs play a central role in the occurrence of liver cirrhosis, the identified miRNAs expression in other hepatic cells and their functions in the pathogenesis of liver cirrhosis deserve further elucidation.

CONCLUSIONS
Based on bioinformatics analysis and cell study verification, we identified 13 differentially expressed miRNAs as potential biomarkers of liver cirrhosis. Moreover, we identified 361 regulatory pairs of miRNA-mRNA and 20 hub genes in liver cirrhosis, most of which were involved in collagen and ECM components, immune response, and cell adhesion. These results would provide novel mechanistic insights into the pathogenesis of liver cirrhosis and identify candidate targets for its treatment.