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

Microarray data

The miRNA expression profile GSE59492 was obtained from the Gene Expression Omnibus database (GEO, http://www.ncbi.nlm.nih.gov/geo/). miRNA expression arrays were performed on liver tissues from five patients with ALD-induced cirrhosis (ALD-CH), 4 patients with HCV-induced cirrhosis (HCV-CH), 5 patients with NASH-induced cirrhosis (NASH-CH) and six patients with noninjured livers (Normal) based on the GPL16384 platform (Affymetrix Multispecies miRNA-3 Array).

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₂ 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.

Immunofluorescence staining

LX2 cells cultured on coverslips at the bottom of 24-well plates were fixed with 4% PFA before permeabilization with 0.1% Triton X-100. After blocking with 10% goat serum, the cells were incubated with anti α-SMA (1:200, Cat# NB600-536, Novus Biologicals, USA) overnight at 4 °C followed by incubation with TRITC-conjugated secondary antibodies (1:200, Cat# ZF-0313, ZSGB-BIO, Beijing, China) at room temperature for 30 min. The cell nuclei were then stained with 4′,6-diamino-2-phenylindole (DAPI, Roche, Switzerland). Coverslips were mounted with anti-fading medium and visualized under a fluorescence microscope (Olympus, Japan).

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.

Identification of DEMs and DEGs

In the GSE59492 profile, 223, 124, and 59 DEMs with P < 0.05 and —fold change— ≥ 1.5 were identified by comparing ALD-CH and Normal, HCV-CH and Normal, and NASH-CH and Normal, respectively (Figs. 2A–2C). The candidate DEMs generated by these comparisons were intersected using a Venn diagram (Fig. 2D). As shown in Fig. 2E, a total of 13 DEMs were identified in cirrhotic livers, of which 6 were 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 7 were 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). Detailed information pertaining to these DEMs was listed in Table 2.

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.

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).

Validation in HSCs by qRT-PCR

To validate the transcriptome data, the miRNA expression levels of the identified DEMs were quantified by qRT-PCR in LX2 cells treated with TGF-β1, which is a classical cellular model for the study of liver fibrosis (Tsuchida & Friedman, 2017). The activation of LX2 cells was confirmed by the overexpression of α-SMA, as proven by immunofluorescence staining (Fig. 8A) and Western blotting (Fig. 8B). After TGF-β1 treatment, hsa-miR-146b-5p, hsa-miR-150-5p and hsa-miR-224-3p expression was upregulated, while hsa-miR-1234-3p, hsa-miR-30b-3p, hsa-miR-3162-3p and hsa-miR-890 expression was downregulated. However, no statistically significant difference was observed in the expression of hsa-miR-3135b, hsa-miR-3195, hsa-miR-4725-3p, hsa-miR-548aj-3p, hsa-miR-548x-3p and hsa-miR-548z in the TGF-β1 group compared with the control group.