Integrated Analysis of Immune-Related circRNA-miRNA-mRNA Regulatory Network in Ischemic Stroke

Background Stroke is the leading cause of death and disability worldwide, with ischemic stroke (IS) being the most prevalent type. Circular RNAs (circRNAs) are involved in the pathological process of IS and are promising biomarkers for the diagnosis of IS. However, studies focusing on circRNAs acting as microRNAs (miRNAs) sponges in regulating mRNA expression are currently scarce. Methods In this study, expression profiles of circRNAs (GSE195442), miRNAs (GSE117064), and mRNAs (GSE58294) from the Gene Expression Omnibus (GEO) database were analyzed. Differentially expressed circRNAs (DEcircRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) were identified by R software. The target miRNAs and target genes were predicted by several bioinformatics methods. Then, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEmRNAs. Subsequently, the protein-protein interaction (PPI) network and the competing endogenous RNA (ceRNA) regulatory network were visualized by Cytoscape software. Finally, we further constructed an immune-related circRNA-miRNA-mRNA regulatory sub-network in IS. Results A total of 35 DEcircRNAs, 141 DEmiRNAs, and 356 DEmRNAs were identified. By comprehensive analysis of bioinformatics methods, we constructed a circRNA-miRNA-mRNA regulatory network, including 15 DEcircRNAs, eight DEmiRNAs, and 39 DEmRNAs. FGF9 was identified as an immune-related hub gene. Immune cell analysis indicated a significantly higher level of neutrophils in IS, and the expression of FGF9 was significantly negatively correlated with the level of neutrophils. Eventually, miR-767-5p was predicted as the upstream molecules of FGF9, and circ_0127785 and circ_0075008 were predicted as the upstream circRNAs of miR-767-5p. Conclusion Our study provides novel insights into the molecular mechanisms governing the progression of IS from the perspective of immune-related ceRNA networks.


INTRODUCTION
Stroke is the primary cause of mortality and disability worldwide (1). Although the stroke mortality rate has declined in recent years, the prevalence has increased due to population growth and aging (2). It is a sudden disruption of cerebrovascular circulation and includes two subtypes: ischemic stroke (IS) and hemorrhagic stroke (3). Ischemic strokes are more common, accounting for ∼87% of all strokes (4). The most effective treatment for IS is intravenous thrombolysis using recombinant tissue plasminogen activator (rTPA). However, the disadvantage of this treatment option is that there is only a 3-4.5 h treatment window (5). Presently, diagnosing IS mainly relies on typical clinical symptoms and brain imaging examinations (6). However, about 50% of IS lack specificity in early imaging diagnosis (7). In addition, the underlying mechanisms of stroke have not been entirely identified, and there are currently no blood biomarkers of IS for clinical use (8). Therefore, it is crucial to explore novel biomarkers of IS, which can contribute to a comprehensive understanding of the etiology and pathophysiology of the disease and aid in early diagnosis and treatment. Previously, miRNA and lncRNA have been proposed as possible stroke biomarkers (9,10). In contrast, researchers have recently begun to investigate the role of circRNA in IS and search for IS-related circRNA biomarkers (11).
Circular RNA (circRNA) is a highly conserved and stable noncoding RNA with a circular structure that does not contain 5′ end caps and 3′ end poly (A) tails (12). CircRNAs have been reported to be extensively involved in the process of IS, including participating in ischemic brain injury, inhibiting apoptosis and neuroinflammation, and protecting the blood-brain barrier (13). In addition, exosomes can carry circRNA across the bloodbrain barrier, allowing circRNA to exist stably in peripheral blood, making circRNA promising as a novel diagnostic and prognostic biomarker for IS (14). A recent study reported that CircOGDH expression was upregulated in the plasma of IS patients, and the CircOGDH expression was positively correlated with penumbra size (15). Li et al. examined circulating blood circRNA expression profiles of acute stroke patients and healthy controls and identified differentially expressed circRNAs in IS patients (16). These previous studies have indicated that peripheral circRNA expression levels are dysregulated in IS patients. Zuo et al. reported that down-regulation of circCDC14A in peripheral blood cells alleviates brain injury in the acute phase of stroke in tMCAO mice (17). Thus, modulating the level of dysregulated circRNAs might reverse the pathological process of IS. circRNA can perform the corresponding function as a competitive endogenous RNA (ceRNA), which means circRNA can bind competitively with miRNA like a sponge through base complementarity, thus affecting the target gene and regulating the level of encoded transcripts (18). A recent study reported that hsa_circ_0045932 (circUSP36) attenuates ischemic stroke injury through the miR-139-3p/SMAD3/Bcl2 signaling axis (19). Xu et al. (20) revealed that circSKA3 could compete for binding hsa-miR-6796-5p and thereby regulates MMP9 to promote ischemic stroke progression. However, there remain numerous dysregulated circRNAs whose roles in IS have not been reported.
Thus, the potential role of circRNAs acting as ceRNAs involved in the process of IS is currently unclear.
Our study conducted a comprehensive analysis of the IS microarray data from the Gene Expression Omnibus (GEO) database. Eventually, we identified a novel potential circRNA-miRNA-mRNA regulatory network in IS. In addition, we screened out immune-related hub genes in the ceRNA network and further constructed an immune-related ceRNA sub-network. Our study provides a novel perspective on the molecular mechanisms of ceRNA immunoregulation in IS.

Microarray Data
Microarray datasets (GSE195442, GSE117064, and GSE58294) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) (21). The dataset GSE195442 is based on the GPL31275 Agilent-085499_SBC human ceRNA microarray [ProbeName version]. All samples in the GSE195442 were included in our study, containing plasma exosome samples from 10 IS patients and 10 controls. The dataset GSE117064 is based on the GPL21263 3D-Gene Human miRNA V21_1.0.0. It contains 1785 serum samples, out of which 189 male participants were collected from individuals aged >65 years, comprising 84 IS patients and 105 controls, and were included in our study. The dataset GSE58294 is based on the GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. Peripheral blood samples from 23 ischaemic stroke patients (5 h of onset) and 23 controls in GSE58294 were included in our study. In addition, GSE16561, based on GPL6883 Illumina HumanRef-8 v3.0 expression beadchip, containing peripheral blood samples from 39 IS patients and 24 controls, was applied for validation. The ethics committee approval or informed consent was not required in this study, as the data were publicly obtained from the GEO database.

Differentially Expression Analysis
DEcircRNAs, DEmiRNAs, and DEmRNAs were screened by the "limma" package in R software (22). The threshold was set with p value <0.01 and |log 2 fold change| >1.2 to select DEcircRNAs from GSE195442. Next, the GSE117064 was analyzed with adjusted p value <0.05 and |log 2 fold change| >2 and set as the threshold for selecting DEmiRNAs. Finally, DEmRNAs were screened by the cut-off point of adjusted p value <0.05 and |log 2 fold change| >1 in GSE58294.

Gene Ontology and KEGG Enrichment Analysis
The "clusterProfiler" package in R was applied for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (23). p-Values lower than 0.05 were considered statistically significant.

PPI Network Construction and Hub Genes Identification
A PPI network of the DEmRNAs was established by the STRING database (http://string-db.org) (24). Interaction scores over 0.4 were considered statistically significant. Next, the interaction network was downloaded and imported into the Cytoscape software for visualization (25). Then, the top 20 genes ranked by Degree and MCC score were selected respectively by the cytoHubba plugin app (26).

Construction of ceRNA Network in IS
We predicted the target miRNAs interacted with the DEcircRNAs by the Circular RNA Interactome (CircInteractome) database (https://circinteractome.nia.nih.gov/) (27). The predicted miRNAs were intersected with the DEmiRNAs by a Venn diagram via the online tool -Venny 2.1.0 (https://bioinfogp. cnb.csic.es/tools/venny/index.html). Then, the target mRNAs of the overlapping miRNAs were predicted by the miRDB database (28), and the intersection of the predicted mRNAs and DEmRNAs was obtained by Venny 2.1.0. Finally, the ceRNA network was constructed based on the expression of circRNAs, miRNAs, and mRNAs, and the results were visualized using Cytoscape software.

Construction of Immune-Related ceRNA Sub-network
Immune-related genes (IRGs) were obtained from the ImmPort database (29). We intersected the IRGs, top 20 hub genes ranked by Degree, top 20 hub genes ranked by MCC, and target genes in the ceRNA network and identified immune-related hub genes. Then, we constructed an immune-related circRNA-miRNA-mRNA sub-network, which was visualized through the Cytoscape software.
Diagnostic Analysis of mRNA of Immune-Related Sub-ceRNA Network An independent external dataset, GSE16561, including 39 IS patients and 24 controls, was used for validation. The effectiveness of the target genes in the immune-related ceRNA sub-network was performed by the receiver operating characteristic (ROC) curves via the "pROC" package in R (30). In addition, the mRNA expression levels of immune-related hub genes were compared between IS patients and controls by t-test, and a p value < 0.05 was considered statistically significant.

Estimation of the Subtype Distribution of Immune Cells
The relative expression of different immune cell subtypes in 46 samples (23 controls vs. 23 IS patients of 5 h onset) in GSE58294 was assessed with CIBERSORT in R (31). First, the proportion of each kind of immune cell in the 46 samples was shown by histogram. Then, the relative expression of every immune cell subtype was compared between the controls and IS group using a boxplot. Finally, the relationship between the immune-related hub genes and each immune cell was performed using Pearson correlation and visualized by a lollipop graph.

Identification of DEcircRNAs, DEmiRNAs, and DEmRNAs in IS
The flow diagram for the whole study is shown in Figure 1. A total of 35 DEcircRNAs (20 up-regulated and 15 down-regulated) were obtained in the GSE195442 dataset (Figures 2A,B). In addition, 141 DEmiRNAs (99 up-regulated and 42 down-regulated) were identified in the GSE117064 dataset (Figures 2C,D). Furthermore, 356 DEmRNAs (165 up-regulated and 191 down-regulated) were screened in the GSE58294 dataset (Figures 2E,F). Lists of DEcircRNAs, DEmiRNAs, and DEmRNAs are available in Supplementary Tables S1-S3.

Gene Ontology and KEGG Enrichment Analysis
Gene Ontology terms are classified into three main categories: biological processes (BP), cellular components (CC), and molecular functions (MF). The top three terms in the BP contain neutrophil activation, neutrophil degranulation, and neutrophil activation involved in immune response. The top three terms in the CC contain external side of plasma membrane, specific granule, and tertiary granule. The top three terms in MF contains: channel activity, ion channel activity, and glycosaminoglycan binding ( Figure 3A). The results of the KEGG enrichment pathway demonstrate that DEGs are principally involved in cytokine-cytokine receptor interaction, PI3K-Akt signaling pathway, and complement and coagulation cascades ( Figure 3B).

Construction of the circRNA-miRNA-mRNA Regulatory Network
A total of 253 potential target miRNAs of the DEcircRNAs were predicted in the CircInteractome database. Then, nine overlapping miRNAs were obtained from the intersection of predicted miRNAs and DEmiRNAs ( Figure 5A). Next, by using the miRDB database, 3,693 potential target genes for the overlapping miRNAs were identified. Subsequently, 72 intersecting mRNAs were acquired by the intersection of predicted target genes and DEmRNAs ( Figure 5B). Finally, the ceRNA network of IS was constructed based on the negative regulation of the circRNA-miRNA pair and miRNA-mRNA pair, which included 15 circRNAs (Figure 5C), eight miRNAs (Figure 5D), and 39 mRNAs (Figure 5E). The ceRNA network is shown in Figure 5F.

Immune-Related ceRNA Sub-network Construction
To further explore immune-related ceRNA regulatory relationships in IS, we obtained immune-related genes (IRGs) in the ImmPort database and intersected the IRGs with mRNAs in the ceRNA network and the hub genes of DEmRNAs ( Figure 6A). As a result, the overlapping gene FGF9 was identified as the immune-related hub gene, and the immune-related ceRNA subnetwork in IS was subsequently constructed. Finally, hsa_circ_0127785 and hsa_circ_0075008 (down-regulated in IS), hsa-miR-767-5p (up-regulated in IS), and FGF9 (down-regulated in IS) formed an immune-related ceRNA regulatory network (Figure 6B). In the independent dataset GSE16561, the expression of FGF9 in IS patients and controls was validated (Figure 7A), and FGF9 expression was significantly associated with a diagnosis of IS (AUC = 0.844) by ROC curve analysis ( Figure 7B).

Relative Immune Cells Expression
We analyzed the relative expression of 22 immune cells in each of the included samples of GSE58294 through CIBERSORT, and the result showed that 17 kinds of immune cells were expressed in the included samples. The distribution of immune cell subtypes in each sample is shown in Figure 8A. Neutrophils recorded the highest expression in all samples, and the expression was significantly higher in IS patients than in controls ( Figure 8B). In addition, the expression of the immune-related hub gene, FGF9, was analyzed in correlation with the expression of each immune cell subtype. As a result, the expression of FGF9 was significantly negatively correlated with the level of neutrophils (p < 0.001; Figures 8C,D).

DISCUSSION
Stroke carries a severe burden on society and the economy (32). IS is the predominant type of stroke. Due to the narrow window for thrombolytic therapy and the frequent lack of specificity of imaging characteristics, it is essential to identify blood biomarkers for rapid diagnosis (33). In addition, identifying potentially relevant genes in stroke and screening for circRNA-miRNA-mRNA regulatory interactions contribute to further understanding of the pathogenesis of IS and provide new insights into the prevention, diagnosis, and treatment of IS.
Numerous studies revealed that abnormal expression of circRNAs is involved in various diseases (34). Recent studies have shown that circRNAs are involved in post-transcriptional regulation in IS (35). Yang et al. (36) revealed that circSCMH1  targeted delivered to the brain mediated by extracellular vesicle improved functional recovery after stroke in mice and monkeys. A clinical trial discovered significant upregulation of circPDS5B and circCDC14A in IS patients (37). Therefore, circRNA is a promising novel target for the diagnosis and prognosis of IS. CircRNA can sponge miRNA, thus reducing miRNA regulation on downstream target genes (38). However, the mechanism of ceRNA in IS is not fully understood, and an increasing number of studies are focusing on the role of circRNA-miRNA-mRNA in IS (39). There is plenty of space to extend further studies on the role of circRNA as ceRNA in IS.
In this study, we identified 35 DEcircRNAs, 141 DEmiRNAs, and 356 DEmRNAs in the expression profiles from the GEO database. Then, we screened circRNA-miRNA interaction pairs through the CircIteractome database and miRNA-mRNA interaction pairs through the miRDB database. By intersecting the predicted miRNAs and DEmiRNAs, and intersecting the predicted mRNAs and DEmRNAs, were finally identified a potential circRNA-miRNA-mRNA regulatory network in IS, which contained 15circRNAs (eight up-regulated and seven downregulated), eight miRNAs (five up-regulated and three  To understand the function of DEmRNAs in IS, we performed the Gene Ontology and KEGG pathway enrichment analysis. As a result, neutrophil activation was significantly enriched in the BP module of GO, and cytokine-cytokine receptor interaction was significantly enriched by KEGG pathway enrichment. Enrichment results confirm the reliability of DEmRNAs, as numerous studies have suggested that the inflammatory response has an essential role in the various phases after the onset of IS, which affects the prognosis and treatment of stroke (40). Thus we further downloaded the immune-related genes (IRGs) from the ImmPort database to screen for the immune-related hub genes in IS. As a result, FGF9, as an intersection of IRGs, PPI hub genes, and target genes in the ceRNA network, was identified as a target gene for constructing immune-related ceRNA sub-network. The sub-network contains hsa_circ_0127785/hsa_circ_0075008 (down-regulated expression in IS), hsa-miR-767-5p (up-regulated expression in IS), and FGF9 (down-regulated expression in IS). The expression pattern of the target gene FGF9 was validated in the independent dataset GSE16561. As a result, FGF9 remained significantly down-regulated expression in IS patients, and the ROC curve indicated a good diagnostic efficacy (AUC = 0.844).
The CIBERSORT algorithm revealed a significant increase in the level of neutrophils after the onset of IS, consistent with previous studies. Studies indicated that properties of circulating neutrophils significantly increased after IS, and the neutrophil alterations were notably correlated with the IS severity (41). An elevated neutrophil to lymphocyte ratio may signal a poor prognosis for IS (42). However, As one of the earliest immune cells to be recruited into the ischaemic brain, neutrophils are involved in a double-edged role in the pathology of IS (43). In addition, the expression of FGF9 was significantly and negatively correlated with the level of neutrophils by Pearson correlation coefficients, which suggested an essential effect of FGF9 in ischemic stroke.  Interestingly, certain recent studies indicated the potential role of FGF9 in IS. It has been reported that long noncoding RNA SNHG7 can mediate up-regulation of FGF9 to alleviate the oxygen and glucose deprivation-induced neuron cell injury (44). Downregulation of FGF9 mediated by miR-339 promotes hypoxia-induced neuronal apoptosis and impairs cell viability (45). These results suggest that FGF9 appears to be a promising target for the diagnosis and treatment of IS. Multiple studies have reported the involvement of hsa-miR-767-5p in cancer. Feng et al. reported that over-expression of hsa-miR-767-5p promotes tumor progression in multiple myeloma model mice, and this effect can be suppressed by CircRNA circ_0000190 (46). Meng et al. suggested that LINC-PINT downregulates miR-767-5p to induce TET2 expression, thus suppressing the aggressiveness of thyroid cancer. However, the role of miR -767-5p in IS has been rarely reported. Herein, we identified a novel immune-related circRNA-miRNA-mRNA regulatory network. We hypothesized that circ_0127785/ circ_0075008/miR-767-5p/FGF9 is involved in the process of ischemic stroke and brings novel insight into the diagnosis and treatment of IS. However, the findings were obtained based on bioinformatic analysis. The regulatory relationships in the ceRNA network remain experimental confirmation.

CONCLUSION
We identified a potential circRNA-miRNA-mRNA regulatory network in IS, including 15 DEcircRNAs, eight DEmiRNAs, and 39 DEmRNAs. In addition, circ_0127785/ circ_0075008/miR-767-5p/FGF9 was identified as an immune-related regulatory sub-network of IS. Our study provides new insights into the pathological processes of IS mediated by circRNA, and the current findings require validation in the future.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
SC designed the study and wrote the manuscript. YZ assisted in analyzing the data and revising the manuscript. MC and WO critically read and edited the manuscript. All authors contributed to the article and approved the submitted version.