Differential Expression Profiles of mRNA and Noncoding RNA and Analysis of Competitive Endogenous RNA Regulatory Networks in Nonalcoholic Steatohepatitis

Nonalcoholic steatohepatitis (NASH) is a liver disease caused by multiple factors, and there is no approved pharmacotherapy. The pathogenesis of NASH remains underexplored. In this study, differentially expressed circular RNAs (circRNAs) were obtained by analyzing NASH-related circRNA datasets, and then, corresponding target microRNAs (miRNAs) and messenger RNAs (mRNAs) were predicted to construct a circRNA–miRNA–mRNA regulatory network. On this basis, a total of 38 circRNAs, 7 miRNAs, and 10 mRNAs were screened out. The present study reveals novel circRNA biomarkers of NASH and reports a potential competing endogenous RNA (ceRNA) regulatory network that might provide insights for further investigation into the underlying pathogenesis of NASH.


Introduction
Nonalcoholic fatty liver disease (NAFLD) is the most common liver disease worldwide and includes nonalcoholic fatty liver, nonalcoholic steatohepatitis (NASH), and cirrhosis. Nonalcoholic fatty liver (NAFL) is characterized by simple steatosis, whereas NASH is typically characterized by the presence of lobular inflammation and ballooning with or without perisinusoidal fibrosis in addition to steatosis [1]. NAFL is the nonprogressive form of NAFLD, while NASH is the progressive form of NAFLD and may advance to cirrhosis and hepatocellular carcinoma (HCC), which is the leading cause of end-stage liver disease or liver transplantation [2]. The prevalence of NASH has been gradually increasing worldwide in recent years, and worryingly, the liver-specific mortality rate for NASH is high [3]. However, the pathogenesis of NASH has not been fully elucidated.
In recent years, various noncoding RNAs (ncRNAs) acting as competing endogenous RNAs (ceRNAs) have become a major research hotspot for various diseases. MiRNAs and circRNAs are different kinds of ncRNAs [4]. Multiple lines of evidence indicate that other RNAs with miRNA target sites, such as circRNAs, can compete with mRNAs to bind miRNAs [5]. CircRNAs have become a focus of life science and medical research and have been identified as key regulators of many diseases. Studies have shown that circRNAs can act as ceRNAs or miRNA sponges by interacting with miR-NAs to sequester these molecules and reduce their regulatory effect on target mRNAs [6]. The circRNA-miRNA-mRNA axis has also been shown to be involved in a variety of cellular events, including apoptosis, vascularization, and metastasis. Studies have shown that the expression profile of circRNA can be a candidate for NASH diagnosis, and the circRNA-miRNA-mRNA pathway may provide clues for studying the pathogenesis of NASH [7,8]. The exploration of circRNA expression patterns and the circRNA-miRNA-mRNA network in the pathogenesis of NAFLD has gradually been carried out. As the pathogenesis of NASH has not been fully elucidated, it appears to be multifactorial. Moreover, the clinical options for NASH are very limited, and many of the drugs in development have failed in both phase 2 and 3 clinical trials. Therefore, research on NASH still faces great challenges. The role of circRNAs in NASH is a new research field. There are few reports on cir-cRNAs in NASH, so further research is needed. Exploring ncRNAs in NASH may provide useful clues to the pathogenesis of NASH.
Therefore, in this study, bioinformatics methods were used to analyze differentially expressed genes (DEGs) associated with NASH, and then, a ceRNA regulatory network involving circRNA, miRNA, and mRNA was constructed to explore some new circRNAs that might be used as ceR-NAs to regulate gene expression in NASH.

Data Collection and Differentially
Expressed circRNA (DEC) Identification. The Gene Expression Omnibus (GEO) is a public functional genomics data repository that supports MIAME-compliant data submissions. This database accepts data based on arrays and sequences. Tools are provided to help users query, locate, review, and download research and gene expression profiles [9]. We searched the dataset of NASH in the GEO database, and a series of related microarray datasets that provide circRNA expression profile data in NASH were acquired. We found raw microarray data for the circRNA expression profile GSE134146 and related GPL microarray gene annotation files [10]. DEC data were obtained by analyzing 4 cases of NASH and 4 controls included in the raw files of the GSE134146 dataset. All raw expression data were normalized by log2 transformation. Then, the online analysis tool GEO2R was used to analyze the differences in microarray data, and the DECs of the microarray dataset were determined with P < 0:05, Log2fold change ðFCÞ > 1 or Log2-fold change ðFCÞ < −1 as the criteria.

Prediction of miRNAs.
CircInteractome computationally identifies potential binding sites for RNA-binding proteins within circRNAs [11]. It maps RNA-binding protein (RBP) and miRNA-responsive element (MRE) sites on human cir-cRNAs by searching some public databases of circRNAs, miRNAs, and RBPs. It uses the TargetScan prediction tool to predict miRNAs that may target circRNAs. miRNet is an easy-to-use web-based tool designed to create, customize, visually explore, and functionally interpret miRNA target interaction networks. It can be integrated into a powerful network visualization system by integrating multiple highquality miRNA target data sources and advanced statistical methods [12]. The relevant target miRNAs of these selected DECs were predicted using two network tools, miRNet and CircInteractome. Overlapping miRNAs for both algorithms were predicted as potential target miRNAs for DECs.
The expression dataset GSE33857 of NASH-related miRNAs was retrieved from the GEO database and includes information for 7 NASH cases and 12 controls. The GSE33857 chip was analyzed by GEO2R, and the differentially expressed miRNAs that overlapped with the predicted targets were included in the next analysis.
2.3. Forecasting of miRNA-Targeted Genes. miRWalk is a comprehensive miRNA target gene database that includes the miRNA target gene information of humans, mice, rats, dogs, cows, and other species. It not only includes the fulllength gene sequence record of the complete miRNAbinding site but is also compatible with the 12 existing miRNA target prediction programs (DIANA-microTv4.0, DIANA-microT-CDS, miRanda-rel2010, mirBridge, miRDB4.0, miRmap, miRNAMap, doRiNA, i.e., PicTar2, PITA RNA22v2, RNAhybrid2.1 and Targetscan6.2). This database can be used to predict the associated combined information set. The two databases miRWalk and miRNet were used to predict target mRNAs of differentially expressed miRNAs. Then, the overlapping DEGs were selected. The GSE24807 dataset is related to the gene expression profiles of NASH, and the raw data were also analyzed with GEO2R [13,14]. Only the DEGs obtained by intersecting the genes of the dataset with the predicted target genes were included in the ceRNA network.

Functional Enrichment Analysis of Overlapping Genes
and Establishment of the Protein-Protein Interaction (PPI) Network. The Database for Annotation, Visualization, and Integrated Discovery (DAVID) provides a comprehensive set of functional annotation tools for investigators to understand the biological meaning behind large lists of genes [15]. It was used to perform Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The Search Tool for the Retrieval of Interacting Genes database (STRING) provides credible information on interactions between proteins and supplies detailed annotation [16].
2.5. Construction of circRNA-miRNA-mRNA Network. To reveal the relationships among circRNAs, miRNAs, and mRNAs, a circRNA-miRNA-mRNA network was constructed by combining circRNA-miRNA interactions with miRNA-mRNA interactions using Cytoscape.

Identification of DECs in NASH.
To construct the interaction network between circRNAs and miRNAs in NASH, DECs should be determined first. A microarray dataset GSE134146 was obtained from the GEO database, which includes 4 NASH cases and 4 controls. The gene chip was from the Agilent-074301 Arraystar Human CircRNA microarray V2 platform. The online analysis tool GEO2R was used to analyze a series of differentiated circRNAs, such as boxplots (Figure 1(a)). P < 0:05 and jlogFCj > 1 were used as the screening criteria in the GSE134146 dataset. There were 192 DECs, including 96 upregulated circRNAs and 96 downregulated circRNAs, between NASH patients and controls ( Figure 1(b)). The expression differences of the top 50 circRNAs with the most significant differences in 4 NASH tissues and 4 control tissues are shown in Figure 2 3.2. Identification of 54 circRNA-miRNA Interactions. Growing evidence indicates that circRNAs regulate gene expression via miRNA sponges. Therefore, some miRNAs related to the DECs we obtained were predicted based on this ceRNA theory. We collected and explored their potential miRNAs through the CircInteractome and miRNet online databases. Among them, 603 miRNAs were found in CircInteractome, and 640 miRNAs were found in miRNet. To ensure accuracy, we used the intersection of the two to obtain 58 overlapping predicted miRNAs. The gene chip dataset GSE33857 from the GEO database was used to verify the predicted miRNAs, and 8 miRNAs that interacted with circRNAs were obtained. According to the related regulatory relationship between circRNA-miRNAs, only 39 circRNAs (25 upregulated, 14 downregulated) and 8 miRNAs (3 upregulated, 5 downregulated) were included in the ceRNA network study. A total of 54 circRNA-miRNA interactions were identified (Table 1).

3.3.
Analysis of miRNA-mRNA Interactions. We obtained 8 miRNAs associated with circRNAs. To explore the functions of these miRNAs in NASH, we used two databases, miR-Walk and miRNet, to predict miRNA-related target genes. A total of 2738 folded predicted target genes were found in both databases. The GSE24807 dataset from the GEO database was used to verify the DEGs. A total of 3245 DEGs were obtained from the dataset. In addition, as shown in Figures 3, 448 overlapping genes were identified by intersect-ing miRNA target genes with DEGs in GEO. Based on the regulatory relationship between miRNAs and mRNAs, 291 genes were included in the list of ceRNAs. A circRNA-miRNA-mRNA regulatory network was constructed by using Cytoscape software (Table 1).

Functional Enrichment Analyses and PPI Network
Construction. Terms related to the DEGs were divided into three functional groups, including biological processes (BP), molecular functions (MF), and cell compositions (CC), using DAVID. The values of the individual components in the GO analysis are shown in Figure 4. In the BP category, 291 DEGs were mainly involved in negative regulation of transcription from RNA polymerase II promoter, positive regulation of transcription from RNA polymerase II promoter, postembryonic development, wound healing, positive regulation of protein kinase B signaling, vascular endothelial growth factor receptor signaling pathway, regulation of defense response to virus by virus, positive regulation of cell proliferation, cell migration, cell motility, and other processes. In the MF category, the genes were mainly enriched in protein binding, growth factor binding, steroid hormone receptor activity, transcription factor binding, ATP binding, transcriptional activator activity, transcription factor activity, sequence-specific DNA binding, DNA binding, double-stranded DNA binding, 1-phosphatidylinositol-3-kinase activity, etc. In the CC category, the genes were mainly enriched in the nucleus, nucleoplasm, cytosol, Golgi apparatus, extracellular exosome, cytoplasm, plasma membrane, lamellipodium, cyclin-dependent protein kinase holoenzyme complex, chromatin, etc. KEGG signaling pathway showed that genes were mainly enriched in cellular senescence, human T cell leukemia virus 1 infection, PI3K-Akt signaling pathway, proteoglycans in cancer, adherens junction, p53 signaling pathway, osteoclast differentiation, pancreatic 3 Gastroenterology Research and Practice cancer, JAK-STAT signaling pathway, Wnt signaling pathway, ErbB signaling pathway, colorectal cancer, lipid and atherosclerosis, pathogenic Escherichia coli infection, oocyte meiosis, cGMP-PKG signaling pathway, T cell receptor signaling pathway, pathways in cancer, cholinergic synapse, axon guidance, MAPK signaling pathway, VEGF signaling pathway, Cushing syndrome, natural killer cell-mediated cytotoxicity, purine metabolism, TGF-beta signaling pathway, focal adhesion, prostate cancer, FoxO signaling pathway, viral carcinogenesis, etc. The results in the KEGG analysis are shown in Figure 5.
3.5. Construction of a circRNA-miRNA-mRNA Network. To further explore the effect of the circRNA-miRNA regulatory network on the expression levels of NASH genes, a PPI network was constructed, and 558 pairs of genes with interactions were found through the STRING database. The PPI network was imported into Cytoscape, and the cytoHubba plug-in was used to further screen hub genes according to the maximal clique centrality (MCC) algorithm. Then, a subnetwork with 10 nodes and 31 edges was selected, which revealed the critical roles of the ten genes (KDR, FYN, RAC1, MAPK1, ERBB2, CDKN1A, HSPA4, SMAD2, MCL1, and ESR1) in NASH ( Figure 6). According to the negative regulatory relationship between ceRNAs, a total of 10 genes and miRNAs were included in the network. After this, a network about the association among these circRNA, miR-NAs, and hub genes was built (Figure 7

Discussion
We successfully constructed a circRNA-related ceRNA regulatory network by integrating and analyzing the expression differences of NASH-related circRNAs, miRNAs, and mRNAs in the GSE134146, GSE33857, and GSE24807 datasets in the GEO database. We found that 39 circRNAs may indirectly regulate 291 mRNAs (or genes) through competitive binding with 8 miRNAs. Among the regulated genes, the 10 most critical central genes were screened out. Then, a net-work of circRNAs, miRNAs, and hub genes was constructed, which contained 38 differentially expressed circRNAs, 7 miR-NAs, and 10 hub genes. These abnormally expressed ceRNAs in NASH have the potential to be excellent biomarkers. The importance of NASH is self-evident, as it may promote the occurrence and development of HCC. NAFLD is a pathological manifestation of metabolic syndrome in the liver. Specifically, it is a form of hepatic steatosis caused by accumulation of liver fat and is closely related to metabolic disorders such as obesity, type 2 diabetes, insulin resistance, hypertension, and hyperlipidemia. NASH is the progressive form of NAFLD. Around 20%-27% of the NAFLD patients develop NASH [17]. NASH is characterized by hepatic steatosis, inflammation, hepatocyte damage, and fibrosis, with inflammation playing a key role in its progression. Liver inflammation is a critical factor in the transition from NAFLD to NASH. Therefore, inflammation is a key pathophysiological mechanism of NASH and a target for therapeutic intervention. The KEGG pathway enrichment results in the current study showed that 291 DEGs were mainly involved in the PI3K-Akt signaling pathway, JAK-STAT signaling pathway, Wnt signaling pathway, cGMP-PKG signaling pathway, T cell receptor signaling pathway,  After multiple screenings, a total of 10 hub genes related to NASH (KDR, FYN, RAC1, MAPK1, ERBB2, CDKN1A, HSPA4, SMAD2, MCL1, and ESR1) in the circRNA-miRNA-mRNA network were identified. Some of them have been linked to liver-related diseases. For example, Zheng et al. identified CDKN1A as a potential key regulator of NASH via dynamic network analysis and dynamic gene coexpression module analysis [18]. Furthermore, studies have shown that the circRNA MAN2B2 promotes the proliferation of hepatoma cells through the miRNA-217/MAPK1 axis [19], which indirectly supports our results. Other studies have shown that HSPA4 is significantly correlated with the prognosis and immune regulation of HCC. Therefore, HSPA4 might be a potential diagnostic and prognostic biomarker as well as a therapeutic target for HCC [20]. In previous studies, these genes, including both MAPK1 and HSPA4, were not reported to be related to NASH but other liver-related diseases. By analyzing the biological processes of these DEGs, we found that these genes may also play an important role in the pathogenesis of NASH.
The importance of ceRNAs in various diseases is emerging, and some ceRNAs have been found to be associated with NASH. There is also evidence that indicates the important regulatory role of miRNAs in NASH [21,22]. Potential targets of differentially expressed miRNAs were known to play a role in lipid metabolism, cell growth and differentiation, apoptosis, and inflammation. For example, overexpression of miR-142-5p inhibits the progression of NASH by targeting thymic stromal lymphopoietin and inhibiting the   JAK-STAT signaling pathway. Thus, miR-142-5p might be a novel latent target for NASH therapy [23]. This is consistent with our findings. In addition, miRNA-223 ameliorates NASH by targeting multiple inflammatory genes in hepatocytes [24]. The role for miR-296 is to regulate lipoapoptosis by targeting p53 upregulated modulator of apoptosis. Hepatocyte lipoapoptosis is a key mediator of liver injury in NASH [25], which makes our data more convincing. Emerging stud-ies seem to establish miRNAs as excellent noninvasive tools for the early diagnosis and treatment of various stages of liver diseases [26]. Recent studies suggest that circRNA may be involved in the pathogenesis of NASH [27]. For instance, steatohepatitis-associated circRNA ATP5B regulator, a mitochondria-located circRNA, was demonstrated to play an important role in alleviating NASH by reducing mROS output [10]. In addition, antagonizing the circRNA_    Figure 6: Hub genes in the NASH-associated PPI network. The 10 hub genes were identified from the PPI network using cytoHubba. The line indicates an interaction between two genes. 11 Gastroenterology Research and Practice 002581-miR-122-CPEB1 axis could alleviate NASH by restoring the PTEN-AMPK-mTOR pathway [8]. However, to date, there is no authoritative agency-approved therapeutic drug on the market. The complex pathogenesis, disease heterogeneity, diagnostic barriers, and selection of treatment endpoints also bring great challenges to NASH research. Therefore, research on NASH still has a long way to go.
Certain potential limitations existed in our study. Further evidence from both in vivo and in vitro experiments is needed for verification. Further study on the physiopathologic mechanism of NASH is being performed on the basis of the current bioinformatics analysis.
In this study, bioinformatics methods were used to integrate NASH and normal liver tissue gene chips to screen out DECs and then search for corresponding miRNAs and competing mRNAs to provide a reference for further research on the pathogenesis of NASH. These circRNAs, miRNAs, and mRNAs were found to be abnormally expressed in NASH, and they have the potential to be potential biomarkers for NASH screening. They also have the potential to enter routine clinical practice and be used as predictive markers of the response to NASH-targeted therapies.

Conclusions
In this study, we constructed a NASH-related ceRNA network by integrating and analyzing the expression differences of NASH-related circRNAs, miRNAs, and mRNAs in the GEO database. These differentially expressed ceRNAs have the potential to be biomarkers for NASH screening and may provide valuable clues for further research on the pathogenesis of NASH.

Data Availability
All data generated or analyzed in this study are contained in the GEO database and this published article.

Conflicts of Interest
The authors declare no conflicts of interest.