Network analysis for the identification of hub genes and related molecules as potential biomarkers associated with the differentiation of bone marrow-derived stem cells into hepatocytes

The incidence of liver diseases has been increasing steadily. However, it has some shortcomings, such as high cost and organ donor scarcity. The application of stem cell research has brought new ideas for the treatment of liver diseases. Therefore, it is particularly important to clarify the molecular and regulatory mechanisms of differentiation of bone marrow-derived stem cells (BMSCs) into liver cells. Herein, we screened differentially expressed genes between hepatocytes and untreated BMSCs to identify the genes responsible for the differentiation of BMSCs into hepatocytes. GSE30419 gene microarray data of BMSCs and GSE72088 gene microarray data of primary hepatocytes were obtained from the Gene Expression Omnibus database. Transcriptome Analysis Console software showed that 1896 genes were upregulated and 2506 were downregulated in hepatocytes as compared with BMSCs. Hub genes were analyzed using the STRING and Cytoscape v 3.8.2, revealing that twenty-four hub genes, play a pivotal role in the differentiation of BMSCs into hepatocytes. The expression of the hub genes in the BMSCs and hepatocytes was verified by reverse transcription-quantitative PCR (RT-qPCR). Next, the target miRNAs of hub genes were predicted, and then the lncRNAs regulating miRNAs was discovered, thus forming the lncRNA-miRNA-mRNA interaction chain. The results indicate that the lncRNA-miRNA-mRNA interaction chain may play an important role in the differentiation of BMSCs into hepatocytes, which provides a new therapeutic target for liver disease treatment.


INTRODUCTION
A liver transplant is a significant way to treat patients with severe liver damage, such as decompensated cirrhosis, liver failure, and advanced liver cancer [1]. However, there is a scarcity of liver donors, and transplantation is associated with immune rejection and other problems [2]. Over the past few decades, in addition to advances in biological treatment research, molecular biology, cell bioengineering, and the stem cell research, stem cell therapy has emerged as an economic and feasible liver disease treatment for the end-stage liver disease [3], particularly decompensated cirrhosis liver failure and advanced liver cancer, and offers an effective strategy with no limit of supply and demand [3]. Stem cell therapy has broad application prospects in liver disease [4]. Bone marrow-derived stem cells (BMSCs) are widely used as adult stem cells, which originate from the mesoderm [5,6]. Several experiments have demonstrated that BMSCs can differentiate into the cells of the mesoderm lineage, such as osteoblasts [7], adipocytes [8], muscle cells [9], neurons and brain cells [10], cardiomyocytes [11], and hepatocytes [12]. Growing evidence suggests that BMSCs can differentiate into hepatocytes, presenting interesting possibilities for cellular therapy of liver diseases. Previous reports have shown that decreased Wnt signaling contributes to the differentiation of BMSCs into hepatocytes [13]. Another study demonstrated the ability of BMSCs to differentiate into liver cells [14]. In addition, Kang et al. reported that rat BMSCs differentiate into hepatocytes [15]. It has been reported that cytokines, such as HGF and bFGF, are key contributing factors in promoting cell differentiation. Under certain conditions, HGF and bFGF can promote BMSCs differentiation into liver cells for the treatment of advanced liver disease [16,17].
Owing to the emergence and development of RNAsequencing technology, a large number of miRNA, lncRNA, and circRNA have been discovered and utilized [18]. LncRNA and miRNA are the two most important types of ncRNA. MiRNAs exhibit posttranscriptional inhibitory effects in animals and plants by pairing with target mRNAs [19]. LncRNA is a class of ncRNAs that does not encode proteins and whose transcripts are longer than 200 nt [20]. LncRNAs are speculated to regulate protein-coding genes in several ways. Studies have found that lncRNA regulates miRNA in three ways: (1) as a precursor or host of miRNAs; (2) lncRNAs and miRNAs compete for mRNA binding; and (3) LncRNA acts as a molecular sponge by absorbing miRNA, thereby regulating gene transcription and expression. Furthermore, miRNA and lncRNA play a major role in various life activities, and in the occurrence and development of liver disease.
Most importantly, these are potential therapeutic targets and diagnostic biomarkers.
The chip data for this study was obtained from a comprehensive gene expression database, and the differential genes were analyzed using bioinformatics software. Different genes were screened to obtain the hub genes that control the differentiation of BMSCs into hepatocytes and liver development. These hub genes include Cat, Cyp2e1, Pah, Ugt2a3, Acss2, Aldh6a1, Hmgcs2, H6pd, Aldh1a7, Hmgcl, Ugt1a1, Arg1, Otc, Baat, Slco1b2, Onecut1, Hhex, Proc, Cdk4, Il6, Fn1, Erbb2, Ccnd1, and Bmp4. In summary, this study provides potential therapeutic targets for the treatment of liver diseases. Figure 1A shows differences in the data, indicating unstandardized data, while Figure 1B shows roughly well-standardized data. A total of 4402 differential genes were detected in hepatocytes, with 1896 upregulated and 2506 downregulated genes compared with BMSCs.

Principal-component analysis (PCA)
The principal-component analysis (PCA) was utilized to obtain information about the overall composition of the analyzed complex microarray datasets. The three samples in the hepatocyte sample group were closely distributed, indicating their high mutual similarity. The BMSCs group exhibited high repeatability. However, the dispersion between the hepatocyte sample and BMSCs groups was very high and well distinguishable, with distinct differences between them, as shown in Figure 2.

Heat map and volcano plot analysis of DEGs
The abscissa represents the sample data, which can be divided into two categories: BMSCs and hepatocytes. In AGING  addition, the expression of genes in the BMSCs and hepatocytes is shown in Figure 3A.
The log2-fold change difference and the negative logarithm of p values between the volcano map samples of DEGs in the BMSC and hepatocyte samples are indicated on the X and Y axes, respectively, each point representing a single gene with detectable expression in both samples. The down-regulated and up-regulated genes were indicated by blue and red, respectively, and insignificant genes were indicated by gray dots. As compared to BMSCs, 1896 and 2506 genes were upregulated and downregulated, respectively, in hepatocytes ( Figure 3B).

GO analysis of DEGs
All DEGs were analyzed by DAVID software (https://david.ncifcrf.gov/tools.jsp). GO enrichment analysis showed that the following GO terms were included for upregulated genes: oxidation-reduction process (GO:0055114), lipid metabolic process (GO:0006629), metabolic process (GO:0008152), and liver development(GO:0001889) ( Figure 4A To screen out the hub genes during the differentiation of BMSCs into hepatocytes, up-regulated genes in hepatocyte oxidation-reduction process, hepatic metabolism process and liver development, and downregulated genes in positive regulation of cell proliferation were selected because they are closely related to hepatocyte differentiation. The PPI network was constructed with four GO genes ( Table 1) and STRING database (https://www.string-db.org/). The DEGs contained in the four GO terms were respectively uploaded to STRING to construct a PPI network. The data exported from the STRING was screened for the hub genes. Two genes were valuable in the oxidationreduction process of hepatocytes: Cat and Cyp2e1. Seven genes were shown to be hub genes within the metabolic process: Pah, Ugt2a3, Acss2, Aldh6a1, Hmgcs2, H6pd, and Aldh1a7. Nine genes play a critical role in liver development: namely Hmgcl, Ugt1a1, Arg1, Otc, Baat, Slco1b2, Onecut1, Hhex, and Proc; and six genes in suppressing liver differentiation: Cdk4, Il6, Fn1, Erbb2, Ccnd1 and Bmp4 (Table 2). Cytoscape v.3.8.2 software act as a pathway for PPI network mapping and analysis. Degree and betweenness centrality (BC) value represents the numbers of interactions of a particular protein and represents nodes passing the node in the shortest distance, respectively. In the final network shown in Figure 5, genes with larger nodes and darker colors are more essential in the PPI network.

Hub genes detected in BMSCs and hepatocytes by RT-qPCR
To investigate key molecular targets regulating differentiation of BMSCs into hepatocytes, four GO terms (Table 1) related to liver development and differentiation were selected. Twenty-four DEGs were identified as hub genes in the differentiation of BMSCs into hepatocytes. Then, twenty-four hub genes were verified by RT-qPCR. Twenty-three hub genes, namely Cat, Cyp2e1, Pah, Ugt2a3, Acss2, Aldh6a1, Hmgcs2, H6pd, Aldh1a7, Hmgcl, Ugt1a1, Arg1, Otc, Baat, Slco1b2, Onecut1, Hhex, Proc, Cdk4, Fn1, Erbb2, Ccnd1 and Bmp4 were consistent with RNA sequencing results, while IL6 gene was contrary to RNA sequencing results. All the 24 hub genes were differentially expressed. Together, these results suggested that twentyfour hub genes are important to regulate differentiation of BMSCs into hepatocytes ( Figure 6).

miRNA prediction
We predicted the miRNAs corresponding to the twentyfour hub genes in miRWalk (http://zmf.umm.uniheidelberg.de/apps/zmf/mirwalk/predictedmirnagene .html). The screening criteria included the miRNA binding site should be in the 3'-UTR, the minimum nucleotide length of miRNA should be 7-mer, and P <0.05. Five prediction programs, miRanda, miRDB, miRWalk, RNA22, and TargetScan, were selected, and the intersection of their predictions was the final prediction result [21]. The genes and MiRNAs interaction network are shown in Figure 7. MiRNAs targeting at least two genes were counted (Table 3).

DISCUSSION
Liver disease is a serious hazard to human health. Ideal clinical intervention drugs and methods have not been identified yet in any part of the world. The liver transplantation is a conventional and effective   intervention method. However, transplantation cannot meet the clinical needs due to the shortage of highquality liver cells, allograft rejection, and other problems [22]. Therefore, there is an urgent need to develop new methods for liver disease intervention. Currently, stem cell-based cell replacement therapy has attracted worldwide attention. BMSC transplantation provides a new way to intervene liver disease. Several studies have shown that transplanted BMSCs can differentiate into hepatocytes to replace the function of damaged hepatocytes and tissues in liver due to their directional differentiation ability, promoting the recovery of liver injury. However, the mechanism of differentiation of BMSCs into hepatocytes remains unclear. The study of key genes and downstream regulatory mechanisms associated with the differentiation of BMSCs into liver cells has a profound clinical impact, which can provide potential targets for BMSC-based treatment of liver failure.
We uploaded chip data of BMSCs and hepatocytes to Transcriptome Analysis Console software for DEGs Analysis. Total 4402 DEGs were detected in hepatocytes, with 1896 upregulated and 2506 downregulated genes compared with BMSCs. GO analysis was performed on all the DEGs, and the results were as follows: among the upregulated genes, 11.44% were related to oxidation-reduction process, 7.16% to liver metabolism process, and 1.92% were to liver development process. Among the downregulated genes, 5.04% were related to cell proliferation. Oxidation-reduction process, metabolism process, liver development process, and cell proliferation are closely related to differentiation of hepatocytes. The STRING search tool was used to build PPI network, and the hub genes were identified by analyzing the network. Among the hub genes, Cat and Cyp2e1 are related to the oxidation-reduction process of hepatocytes; Cat and FOXO3 are positively correlated with the differentiation of BMSCs into cells of the osteogenic lineage [23]. In the later stage of stem cells differentiation into hepatocytes, Cyp2e1 is expressed [24]. The hub genes Pah, Ugt2a3, Acss2, Aldh6a1, Hmgcs2, H6pd, and Aldh1a7 are related to liver metabolism. The mRNA level of Hmgcs2 increase during the differentiation process of embryonic stem cells to hepatocellular-like cells [25]. The role of Pah, Ugt2a3, Acss2, Aldh6a1, H6pd, or Aldh1a7 in hepatocyte cells is currently unclear. The hub genes Hmgcl, Ugt1a1, Arg1, Otc, Baat, Slco1b2, Onecut1, Hhex, and Proc are related to liver development. Ugt1a1 is a marker of hepatocytes, and its significant expression during the differentiation of human hematopoietic stem cells into hepatocytes indicates that HSCs have successfully differentiated into normal hepatocytes [26]. Arg1 expression is observed during hepatic-like phenotype differentiation of somatic stem cells in vitro [27]. In a previous study, increased Otc expression was reported in 201B7 cells after culture in hepatocyte differentiation initiation medium [28]. Onecut interacts with Lmx1a to promote the differentiation of ventral midbrain neural stem cells into dopamine neurons through the Wnt1-Lmx1a pathway [29]. Hhex regulates the hepatic differentiation of embryonic stem cells [30]. The role of Hmgcl, Baat, Slco1b2, or Proc in hepatocyte differentiation and development is unclear. The hub genes Cdk4, Il6, Fn1, Erbb2, Ccnd1, and Bmp4 regulate cell proliferation networks and are downregulated in hepatocytes compared with untreated BMSCs, suggesting inhibitory effects on liver differentiation. Studies have shown that expression of Cyclin B1 and Cdk4 during the hepatic differentiation of liver epithelial progenitor cells (LEPCs) induced by sodium butyrate may be related to the growth arrest of LEPCs shortly after treatment [31]. MSCs contribute to liver differentiation by activating the IL-6 / gp130-mediated STAT3 signaling pathway [32]. Ccnd1 has been reported to be associated with liver regeneration, and it is speculated that they play a key role in mouse hepatocytes [33]. Ccnd1 silencing suppresses liver cancer stem cells (LCSCs) differentiation [34]. Bmp4 is an important regulator of  cell proliferation and differentiation. Studies have shown that Bmp4 is a key cytokine for the development of mouse embryonic stem cells into hepatocytes [35]. Fn1 and Erbb2 have not been reported in this respect. Furthermore, we verified twenty-four hub genes in the BMSCs and hepatocytes using RT-qPCR.
All the twenty-four hub genes were differentially expressed.
Taken together, through bioinformatics analysis, we identified key genes that regulate the differentiation of BMSCs into hepatocytes and their upstream miRNAs and lncRNAs, providing potential targets for stem-cell replacement therapy for liver diseases. Hence, exploring the key genes and downstream regulation mechanism of the differentiation of BMSCs into hepatocytes is essential for treating hepatic diseases.

Acquisition of primary mouse BMSCs and hepatocytes
In this experiment, two-step perfusion method was used to isolate and acquire mouse primary hepatocytes. First, the liver was flushed with Ca2 + -free HEPES, followed by HEPES perfusion containing collagenase and Ca2 + , and then the cell suspension was cultured in DMEM/Ham's F12 medium. The primary hepatocytes were established in complete DMEM/Ham's F12 media placed in a 95% air/5% CO2, 37° C incubator [41].
The soft tissue on the surface of the mouse femur was carefully peeled off and the femur tissue was rinsed several times with a 1 mL syringe containing α-MEM culture solution, until it was transparent. Four rinsed femur cells were cultured in a medium containing BMSCs growth medium placed in a 95% air/5% CO2, 37° C incubator, and refreshed every three days [42][43][44].

Download of microarray data
Gene Expression Omnibus (GEO) database provides gene expression profiles [45]. We downloaded two gene expression profiles of GSE30419 and GSE72088 from the GEO database datasets. GSE30419 included three untreated Mus musculus BMSC samples (GSM795638, GSM795639, and GSM795640), and dataset GSE72088 included three PBS treated Mus musculus primary hepatocytes samples (GSM1375704, GSM1375705, and GSM1375706).

Screening of DEGs
The genes expression profiles of GSE30419 and GSE72088 were analyzed by Transcriptome Analysis Console software. | Fold Change | >4, P <0.05 and FDR <0.05 were used as the standard for analysis and screening of DEGs.

GO enrichment analysis
To investigate the functions of these gene signatures, we performed GO enrichment analysis. The DAVID website has multiple functional enrichment for DEGs [46], and the screening criteria for GO functional enrichment is usually set as P <0.05 and FDR <0.05 [47]. The results were presented as bar graphs by Graph Pad Prism 6.0.

Discovery of hub genes
The PPI network of differential genes was constructed using the STRING online platform and analyzed and processed by Cytoscape V.3.8.2 software [48]. The means of the degrees were calculated to screen out genes that showed a greater degree than the mean value. The hub genes with BC value greater than 0.05 and degree exceeds the mean were selected.

Verification of hub genes expression using RT-qPCR
Total RNA was extracted from hepatocytes and BMSCs using TRIzol reagent (Sigma, St. Louis, MO, USA). We used the cDNA Synthesis Mix Kit (Innovagene, Hunan, China) to generate the first strand cDNA. Real-time PCR was implemented using SYBR Green qPCR Mix (Innovagene). Using the 2 −ΔΔCq method, we calculated and analyzed the expression situations of hub genes.

Gene-miRNA prediction
We used the miRWalk website to predict target miRNAs located upstream to the genes [49]. Five databases including miRanda, miRDB, miRWalk, RNA22, and TargetScan were selected for prediction. The interaction between genes and miRNAs was visualized by Cytoscape V 3.8.2. MiRNAs targeting at least two genes were screened.

miRNA-lncRNA prediction
Interaction between miRNAs and lncRNAs was identified using StarBase website [50], and Cytoscape V 3.8.2 was used for analysis and mapping. LncRNAs targeting multiple miRNAs are required to be carefully analyzed. AGING

Statistical analysis
Statistical analysis of the data was performed in IBM SPSS Statistics Version 26. The data based on RT-qPCR experiment was expressed as mean ± SD and presented using Graph Pad Prism 6.0. All the genes were tested in three separate experiments. Analysis was done using student's unpaired t-tests. P<0.05 was statistically significant.

AUTHOR CONTRIBUTIONS
YHH, XMH, SJL, and TK contributed to the conception of the study, writing the manuscript, and performing the literature search. YYM, XCL, and SJL performed the data analysis. YHH, MHJ, and HNS performed analysis and the quality assessment of the study. YHH and TK designed the experiments, supervised the research, and prepared the manuscript. All authors read and approved the final manuscript.