Comprehensive Analysis of RNA Expression Profile Identifies Hub miRNA-circRNA Interaction Networks in the Hypoxic Ischemic Encephalopathy

Hypoxic ischemic encephalopathy (HIE) is classified as a sort of serious nervous system syndrome that occurs in the early life period. Noncoding RNAs had been confirmed to have crucial roles in human diseases. So far, there were few systematical and comprehensive studies towards the expression profile of RNAs in the brain after hypoxia ischemia. In this study, 31 differentially expressed microRNAs (miRNAs) with upregulation were identified. In addition, 5512 differentially expressed mRNAs, long noncoding RNAs (lncRNAs), and circular RNAs (circRNAs) were identified in HIE groups. Bioinformatics analysis showed these circRNAs and mRNAs were significantly enriched in regulation of leukocyte activation, response to virus, and neutrophil degranulation. Pathway and its related gene network analysis indicated that HLA − DPA1, HLA − DQA2, HLA − DQB1, and HLA − DRB4 have a more crucial role in HIE. Finally, miRNA-circRNA-mRNA interaction network analysis was also performed to identify hub miRNAs and circRNAs. We found that miR-592 potentially targeting 5 circRNAs, thus affecting 15 mRNA expressions in HIR. hsa_circ_0068397 and hsa_circ_0045698 were identified as hub circRNAs in HIE. Collectively, using RNA-seq, bioinformatics analysis, and circRNA/miRNA interaction prediction, we systematically investigated the differentially expressed RNAs in HIE, which could give a new hint of understanding the pathogenesis of HIE.


Introduction
Hypoxic ischemic encephalopathy (HIE) is classified as a sort of serious nervous system syndrome that occurs in the early life period arising from placental insufficiency or umbilical cord obstruction in the perinatal period [1]. HIE is thought to be the chief neuro-developmental disorder in infants, with an estimated incidence rate of 1 to 6 per 1000 births [2,3]. It has shown that HIE constitutes a primary inducer of incidence and lethality rates in neonates [4]. Afterwards, cell damage in the central nervous system caused by neonatal hypoxia and ischemia is still the common contributor to neurological disorders in adulthood, such as cerebral palsy and mental retardation [1][2][3][4][5]. Despite mounting studies have revealed the morphological, biophysical, and biochemical changes that occur after neonatal hypoxia ischemia (HI) brain damage, the efficacy of neuro-protective treatments is limited to HIE patients [6]. In order to give impetus to neurotherapeutics and ameliorate clinical outcomes in HIE, powerful biomarkers towards brain injury are essential to quickly assess the effect of treatment, offer prognostic details for family members, and guide the needs of later recovery care, thereby contributing to personalized care [6,7].
Currently, miRNAs have attracted the attention of many researchers due to their crucial roles in posttranscriptional regulation in human diseases. There are several studies that have shown most miRNAs are observed in brain tissues. Unique properties of miRNAs as potential markers for a variety of diseases, such as oncology and heat disease [8][9][10]. Additionally, circular RNA (CircRNA) has been confirmed to participate in many biological functions, comprising the promotion of rolling circle translation, the control of parent genes transcription, the assistant of alternatively spliced mRNA formation, and the sponging actor of microRNA (miRNA/miR) [11,12]. It has been shown that changes in the levels of specific circRNAs exhibited an association with carcinoma, ischemia, stroke, neurodegenerative diseases, heart disease, etc. [13,14]. CircRNAs, especially in brain tissue, were not only related to the nervous system's development, differentiation, and biological function but also functioned importantly in the dysfunction and pathology after brain damage [15]. Recently, it has been reported that circRNAs' expression profile in brain tissue was modified under cerebral ischemia, among which were probably involved in the etiopathogenesis of cerebral ischemia. In adult rats, there were some studies [11,16] showing that the expression profiles of cerebral circRNA altered significantly following stroke, which was conducive to stabilizing mRNA expression. Also, several circRNAs had been revealed to be related to HIE. For example, circular RNA cZNF292 silence alleviates neural injury through upregulation of miR-22 in rat model [17]. In addition, Jiang et al. identified a total of 66 circRNAs were differentially expressed in HIBD rats compared with the control group [18]. So far, there were few systematical and comprehensive studies towards the expression profile of RNAs in the brain after hypoxia ischemia.
In our study, we attempted to investigate novel diseaserelated RNAs in HIE-induced rat model to decipher the biological function of RNAs. RNA sequencing (RNA-seq) is a powerful tool for analyzing the changes of transcriptome in cells and tissues [19]. This method can obtain information related to RNA splicing, gene expression, and single nucleotide polymorphism (SNP) change [19]. It can also be utilized as a distinguishing tool for ribosomal RNA, transfer RNAs, microRNAs, and small RNAs. In recent years, RNA-seq has been largely applied to identify the differential expression of genes. Herein, the expression profiles of RNAs were analyzed a by RNA-seq. The possible function of RNAs in HIE was predicted by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Furthermore, miRNAs regulated by differentially expressed RNAs were predicted in TargetScan database. Collectively, the combined use of RNA-seq data and bioinformatics analysis in out literature uncovered the potential roles of these crucial RNAs in pathophysiological processes of HIE.

Materials and Methods
2.1. Sample Collection. The two RNA-seq dataset and the corresponding metadata were acquired from Genes Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih .gov/geo) database (GSE164727 and GSE121178 [20]). Three healthy neonatal rats (control group) and three HIBDinduced rats (experimental group) were distributed in GSE164727 dataset. Three whole blood samples from HIEinduced rats and three normal blood samples from control group were allocated in GSE121178 dataset [20].

Differential Expression
Analysis. The count matrix of gene expression and the corresponding grouping data were obtained from GEO database and converted to count-permillion (CPM) for difference analysis. The limma package in the R software (version 3.3.2, https://www.r-project.org/) was utilized to screen the differentially expressed miRNAs and circRNAs, with the selection criterion of adjusted p value < 0.05 and |log2 fold change ðFCÞ | >1. The DEGs with log2 FC < −1 were DEGs with downregulation, whereas those with log2 FC > 1 were DEGs with upregulation.     3. Enrichment Analysis. Gene Ontology (GO), KEGG pathway analysis, and Reactome pathway analysis were performed. Hypergeometric distribution was applied to calculate the p value. R package "clusterProfiler" (version 3.3.2) was employed to conduct all enrichment analysis [28]. p < 0:05 meant there was greatly significant difference in two or more compared groups.
2.4. Gene Set Enrichment Analysis (GSEA). GSEA is taken to determine whether a set of a priori defined genes display significant statistics differences in two biological states (e.g., phenotypes). GSEA was conducted by the GSEA software (version 1.46.0) [21].
2.6. Statistical Analysis. R (version 3.6.1) was applied for statistical analysis, and Wilcoxon rank sum test was utilized to investigate the difference between the two groups. p < 0:05 meant that significantly statistics difference existed in compared groups. The pheat map package (version 1.0.12) was applied to plot the heat map for the DEGs. The R package "ggbiplot" (version 0.55) was employed to conduct the principal component analysis (PCA) for miRNAs.

Identification and Characteristic
Comparisons of miRNAs, circRNAs, and mRNAs. The flow chart analysis is illustrated in Figure 1. Totally, 627 miRNAs were identified. After excluding miRNAs with low expression (mean expression < 1), 179 miRNAs were retained for subsequent screening and applied for PCA analysis. Figure 2(a) shows that the majority of samples in the same group were classified into the same group, indicating that the miRNA expression between model mice and normal mice was significantly different. 31 differentially expressed miRNAs with upregulation were identified in HIBD/HIE rats group after comparison with the normal group ( Figure 2(b)). We further evaluated the expression value for whole samples with the differentially expressed miRNAs between the two groups, and the most top 12 differentially expressed miRNAs were chosen to construct a boxplot, including rno_mir_1193, rno_mir_133c, rno_mir_186, rno_mir_2985, rno_mir_329, rno_mir_344i, rno_mir_3571, rno_mir_379, rno_mir_592, rno_mir_6327, rno_mir_679, and rno_mir_758 ( Figure 2(c)).
Next, we normalized the expression profile data (Figure 3(a)). The UMAP plot shows that the HIE group and normal group have different expression patterns (Figure 3(b)). A total of 5512 differentially expressed mRNAs, lncRNAs, and circRNAs were identified in the experimental and control groups. Among them, 2363 genes were upregulated compare the HIE group to the normal group, and 3149 genes were downregulated (Figure 3(c)). Then, we normalize the expression values of these 5512 genes and draw the heat map with cluster the samples. The results show that the samples in the identical group can be clustered together well (Figure 3(d)).

Functional Enrichment Analysis for All Differentially
Expressed circRNAs and mRNAs. There are 29 GO terms that were enriched in the differentially expressed circRNAs and mRNAs which include 14 Cellular Component (CC) terms, 4 Biological Process (BP) terms, and 9 Molecular Function (MF) terms (Figure 4(a)). Among them, regulation of leukocyte activation, negative regulation of cell activation and leukocyte activation, and cellular response to hepatocyte growth factor stimulus were the most significantly enriched in BP (Figure 4(b)). The GSEA results showed that 136 GO terms were enriched for the differentially expressed RNAs. The most significantly enriched term in BP is response to virus (Figure 4(c)).
Additionally, there were 8 pathways in KEGG and 25 pathways in Reactome enriched for these genes, respectively (Figures 5(a) and 5(b)). And the pathways for amoebiasis, intestinal immune network for IgA production, allograft rejection, and hematopoietic cell lineage are the most enriched pathways in KEGG (Figures 5(c) and 5(e)). Pathway and its related genes network analysis indicated that the pathways with significant enrichment have a strong correlation ( Figure 5(d)). Multiple pathways were potentially regulated by several hub genes; for example, ICOS, HLA − A, and CD40LG were related to regulate allograft rejection and cell adhesion molecules; IL10 was related to regulate allograft rejection, amoebiasis, and hematopoietic   3.3. miRNA-circRNA-mRNA Interaction Network Analysis. 2432 differentially expressed circRNAs were identified. Among them, 772 have comments in CDC database. A total of 2055 miRNAs was predicted in TargetScan database. Then, we matched these predicted miRNAs with the miR-NAs screened in the previous section, and the results showed that only miR-592 is available in both datasets. Through screening, we found that a total of 15 circRNA (hsa_circ_ 0011142, hsa_circ_0015469, hsa_circ_0025546, hsa_circ_ 0042325, hsa_circ_0045698, hsa_circ_0045932, hsa_circ_ 0058532, hsa_circ_0059400, hsa_circ_0064454, hsa_circ_ 0067953, hsa_circ_0068397, hsa_circ_0077377, hsa_circ_ 0082316, hsa_circ_0084645, hsa_circ_0088562) response elements contain miR-592, and these circRNAs are in different gene regions. And then, we used Cytoscape to draw the network for all miRNA, circRNAs, and mRNAs ( Figure 6(a)). Among these circRNAs, the bioinformatics analysis demonstrated that hsa_circ_0045698 and hsa_circ_ 0068397 were the most significantly circRNAs, which including the most miRNA binding sites. As present in Figure 6(b), the CSCD database was used to draw a circle graph for these miRNAs of which targets are also differentially expressed (Figure 6(b)).

Discussion
HIE is a main contributor to higher incidence and fatality rates in neonates. Accompanied by the advancement of obstetrical and neonatal care, the survival probability of HIE has also heightened [1,3,25]. Nevertheless, the presently available treatment for therapeutic hypothermia can only benefit about 10% of babies. More importantly, it is still elusive towards the mechanisms of HIE. HIE is still considered as a detriment to children's health and life quality, and it is urgent to explore efficacious neurotherapeutic interventions and promising biomarker in the treatment of HIE [26,27].
Until now, there are limited publicly reported studies on the expression of miRNAs in HIE. Shi et al. found that miRNA-21 and HIF-1α expression levels were higher in HIE neonate serum samples using a small cohort of newborns with HIE. Cai et al. showed that overexpression of miR-27a weakened HI-induced neuronal apoptosis via modulating FOXO1 in rat model [28]. Hypoxia-induced miR-152 suppresses cell apoptosis via inhibiting PTEN [29]. The present study determined the expression profile of differentially expressed miRNAs using RNA-Seq along with bioinformatics analysis. A total of 627 miRNAs were detected, 179 of which were chosen for PCA analysis. 31 differentially expressed miRNAs with upregulation were identified after comparing the two groups, and the most top 12 differentially expressed miRNAs, including rno-mir-1193, rno-mir-2985, rno-mir-6327, rno-mir-3571, rno-mir-379, rno-mir-344i, rno-mir-133c, rno-mir-329, rno-mir-758, and rno-mir-592. Among these miRNAs, miR-329-5p was found to interact with upregulated mmu-circRNA-015947 in neuron injury condition [18]. Mir-592 was reported to regulate the induction and cell death-promoting activity of p75NTR in neuronal ischemic injury [30]. miR-344 was reported neural-specific expressed during mouse embryonic development, indicated it had a crucial role in brain development [31]. Abnormal expression of mir-379 implies the pathogenesis of spinal cord injury [32]. These previous reports and this study indicated these miRNAs may play a key role in brain injury.
The microarray analysis also consisted of the expression profile of coding genes. We further evaluated the differentially expressed RNAs between HIE and normal groups. A total of 5512 differentially expressed mRNAs, lncRNAs, and circRNAs were identified in HIE and control groups. Among them, 2363 genes were induced, and 3149 genes were suppressed after comparison with the two groups. Typically, hsa_circ_0079200 and Hsa_circ_ 0083669 were significantly upregulated, while hsa_circ_0073814, hsa_circ_ 058702, and hsa_circ_0070224 were greatly downregulated. LNCV6_97168 was upregulated, whereas LNCV6_132869, LNCV6_56111, and LNCV6_34762 were downregulated. These results suggested that these altered lncRNAs and cir-cRNAs in the pathophysiological processes following HIE  could be potential biomarkers or promising therapeutic targets in treating this disease. GO analysis of these circRNAs and mRNAs was significantly enriched in regulation of leukocyte activation and cellular response to hepatocyte growth factor stimulus. The GSEA data revealed that there are 136 enriched GO terms for the differentially expressed RNAs, and the most significantly enriched biological process term is response to virus. KEGG analysis shows that amoebiasis, intestinal immune network for IgA production, allograft rejection, and hematopoietic cell lineage are the most enriched pathways. GSEA Reactome analysis indicated that neutrophil degranulation, signaling by interleukins, GPCR ligand binding, and cellular response to stress were the most significantly enriched terms for these RNAs. Neutrophils are the first type of peripheral inflammatory cells in acute and chronic inflammation. Several lines of evidences have shown that infiltration of neutrophils in the adult cerebral ischemic area within a few hours after injury. In the neonatal HI injury model, the response of neutrophils in the brain is usually weakened because they seem to stay mainly in the blood vessels [33]. The expression of interleukin 1β (IL-1β) in cord blood and peripheral blood of HIE children was significantly upregulated, and the upregulated level was corresponding to HIE grades and adverse outcomes [34]. The release of IL-1β relied on the activation of NLRP3 inflammasome [35]. In disease models such as cerebral hemorrhage model and ischemia-reperfusion injury model, NLRP3 inflammasome was reported to be motivated in astrocytes [36,37]. Mitochondria are the main target of HI injury. The impairment of mitochondrial increased with age, bringing about the disorder of mitochondrial-related molecular pathways. During HI injury, reactive oxygen species (ROS) generated by mitochondria exceeded the antioxidant capacity, resulting in DNA damage, mitochondrial lipid peroxidation, and Ca2+ homeostasis disruption. Therefore, mitochondrial function played as key determinant of survival during ischemic injury [38].
Multiple genes were dramatically upregulated in neonatal HIE, including IL10, ARG1, and IL1R2. Interleukin 10 (IL-10) is a synthetic human cytokine inhibitor with immunomodulatory and anti-inflammatory effects. It is the key to successfully resist the damage of host tissue during the acute phase of immune responses [39]. A large number of studies have verified the changes of IL-10 as well as its important role in HI injury, indicating that the release of IL-10 was increased, and it had a negative correlation with the rate of neuronal apoptosis [46]. For instance, Bai et al. [39] reported that the level of IL-10 was markedly raised after HI and IL-10 were colocalized with HI insult-injured neurons and astrocytes. IL-10 was shown to exhibit a repressive effect on immune responses and a decreased inflammatory reactions along with neuronal injury [40] in cortical tissues [40] and in serum cytokines of the liver [40]. Despite the expression of IL-10 was largely studied, the mechanism of IL-10 involving in HIE ignition and progression was still elusive [41]. Arginases are pivotal modulatory enzymes of inflammation and tissue repair [42]. At present, it is reported that the expression of ARG-1 in the blood of patients with stroke is heightened, and its level is related to infarct size [43], implying the involvement of ARG in the mechanisms of post-hypoxic-ischemia. Till now, there were emerging researches that have showed that ARG-1 was primarily related to the role of an indicator of "pro-repair" microglia, while ARG-2 was mainly related to cerebral blood flow regulation [44]. Hypoxia and ischemia-reperfusion injury are considered to be powerful stimuli for the expression of ARG-1 and ARG-2. IL1R2 is a cytokine receptor of the interleukin-1 (IL-1) receptor family, and it plays as an essential mediator of several cytokines generated by immune and inflammatory responses. Highly expressed IL1R2 was positively correlated with IL1R signals, manifesting that IL1R2 may participated in the interplay between IL1R1 and IL-1. In hypoxia/reoxygenation-(H/R-) treated cardiomyocytes, overexpressing IL1R2 attenuated the promoting proliferation and antiapoptosis effects induced by the overexpression of LncRNA A2M AS1 [45].
It has been reported that circRNAs act as a miRNA sponge and inhibit miRNA activity, leading to the upregulation of miRNA targets [46]. miRNAs had a crucial role in the central nervous system. Some circRNAs may took part in HIBD development through circRNA/miRNA interactions. miRNA target prediction software was utilized to prove that the dysregulated circRNAs that included miRNA binding sites [47]. A total of 2,432 differentially expressed circRNAs were identified, amid which 772 were annotated in CDC database. 2055 miRNAs were predicted in TargetScan database, and we found that only miR-592 is available in both datasets. Moreover, many previous studies had demonstrated that miR-592 was related to multiple nervous disorders [48][49][50]. A previous report has identified regionspecific differences in miR-592 levels in the brain of adult rodent. Additionally, some reports have presented that miR-592 was differentially expressed in multiple neoplasms. One study showed that the altered miR-592 exerted impressive consequences in the nervous system by modulating the expression of p75NTR [30].
Finally, we found the most significantly differentially expressed circRNAs (hsa_circ_0045698 and hsa_circ_ 0068397). Future researches will concentrate on exploring more downstream targets of both two circRNAs, helping to decipher their role in the pathophysiological process of neonatal HIBD. Collectively, using RNA-seq, bioinformatics analysis, and circRNA/miRNA interaction prediction, we systematically investigated the differentially expressed RNAs in HIE model. These discoveries could give a new hint of understanding the pathogenesis of HIE and in depth deciphering the biological functions of these candidate circRNAs/miRNAs.

Data Availability
The data used during the present study are available from the corresponding author upon reasonable request.

Conflicts of Interest
No conflict of interest.

Authors' Contributions
Lin Wei and XiaLi analyzed the data and wrote the manuscript. Lijuan Wang and Yanyan Song designed the study. Hongmei Dong revised the manuscript. All the authors agreed to be accountable for the accuracy and integrity of all aspects of the research. Lin Wei and Xia Li contributed equally to this work. 16 Computational and Mathematical Methods in Medicine