Identifying circRNA-associated-ceRNA networks in retinal neovascularization in mice

Retinal neovascularization is a complication which caused human vision loss severely. It has been shown that circular RNAs (circRNAs) play essential roles in gene regulation. However, circRNA expression profile and the underlying mechanisms in retinal neovascular diseases remain unclear. In the present study, we identified altered circRNAs in the retinas of oxygen-induced retinopathy (OIR) mouse model by microarray profiling. Microarray analysis revealed that 539 circRNAs were significantly altered in OIR retinas compared with controls. Among them, 185 up-regulated and 354 down-regulated circRNAs were identified. The expression levels of 4 altered circRNAs including mmu_circRNA_002573, mmu_circRNA_011180, mmu_circRNA_016108 and mmu_circRNA_22546 were validated by quantitative real-time reverse transcription-polymerase chain reaction (qRT-PCR). Bioinformatic analysis with validated circRNAs such as competing endogenous RNA (ceRNA) regulatory networks with Gene Ontology (GO) enrichment analysis demonstrated that qRT-PCR validated circRNAs were associated with cellular process, cell part and phosphoric ester hydrolase activity. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis demonstrated that MAPK signaling pathway and renin-angiotensin system were related to validated circRNAs, suggesting these pathways may participate in pathological angiogenesis. The results together suggested that circRNAs were aberrantly expressed in OIR retinas and may play potential roles in retinal neovascular diseases.


Introduction
Ischemia-induced retinal neovascularization is a key pathological process of many retinal vascular diseases like diabetic retinopathy and retinopathy of prematurity (ROP). It is also a leading cause of vision loss as well as blindness in developed countries and regions [1]. Studies have shown that a variety of altered genes are involved in the pathogenesis of retinal neovascularization [2,3]. Several molecules, such as vascular endothelial growth factor (VEGF), basic fibroblast growth factor as well as periostin, have been identified as key mediators in pathological angiogenesis [4,5]. In particular, VEGF has been widely studied and targeted in clinical applications where reports showed anti-VEGF treatment was remarkably efficient in patients [6]. However, some patients are not responsive to anti-VEGF therapy [7], suggesting further studies should aim to look for other therapeutic targets and biomarkers.
Circular RNA (circRNA), as a type of endogenous non-coding RNA, is characterized by covalently closed continuous loop structure lacking 5' cap and poly(A) tail at 3' ends [8]. CircRNA is widely expressed in serum and tumor tissues in mammals [9,10]. CircRNA functions as a molecular sponge efficiently targeting miRNA and inhibiting miRNA transcription [11,12]. Through targeting miRNA, circRNA regulates downstream gene expression and may play a crucial role in disease mechanisms [13]. A

Ivyspring
International Publisher few studies have demonstrated that several circRNAs participate in the pathogenesis of diabetic retinopathy [14,15], vascular endothelial cells proliferation and angiogenesis [16], suggesting circRNAs could serve as potential biomarkers for diagnosis of diabetic retinopathy and provide novel therapeutic targets to treat diabetic retinopathy. A study reported that circRNA-MYLK served as competing endogenous RNA (ceRNA) for miR-29a, contributing to epithelial-mesenchymal transition and the bladder cancer progression by activating VEGFA/VEGFR2 pathway [17]. Besides, circRNAs ZNF280C_hsa_circ_ 001211 and SIAE_ hsa_circ_002083 may participate in the key pathways of ROP pathogenesis [18]. However, the involvement of other potential circRNAs and the underlying mechanisms of ROP remain unclear.
We previously demonstrated differentially expressed mRNAs and long non-coding RNAs in the oxygen-induced retinopathy (OIR) mouse model [19], suggesting that non-coding RNAs played crucial and different roles in retinal neovascularization. In this study, we performed microarray analysis aiming to profile expressions of another kind of emerging non-coding RNA, circRNA in retinas of OIR mice. Furthermore, we interrogated the putative functions of the altered circRNAs by bioinformatic analysis.

Animals and statement of ethics
C57BL/6J mice (Hunan SJA Laboratory Animal, Changsha, China) were used in the present study. The animal experiments were conducted based on ARVO Statement for the Use of Animals in Ophthalmic and Vision Research, and were subjected to approval by the Institutional Animal Care and Use Committee of Central South University.

Oxygen-induced retinopathy mouse model
The OIR mouse model was established as described [20]. Pups were exposed to hyperoxia environment (75% oxygen) for 5 days at postnatal day 7 (P7), followed by return to the environment of room air at P12. Newborn pups kept continuously in the room air environment were used as the controls. The retina samples were collected from mice of both groups at P17 for analysis.

RNA isolation
Total RNAs were isolated from retinas (retinas from both eyes of the same mice were pooled as one sample) by using TRIzol reagent (Invitrogen, Carlsbad, USA). The concentrations were assessed by the NanoDrop ND-1000 (Thermo Scientific, Wilmington, DE, USA). The RNA integrity was examined by electrophoresis on a denaturing agarose gel prior to further investigation.

Microarray analysis of circRNAs
Microarray analysis is widely used to identify altered circRNAs in in vivo studies [21]. Three pairs of retinas in each group were used for the analysis. Sample labeling and microarray hybridization were conducted according to the standard protocols Arraystar (Rockville, MD, USA). Briefly, total RNAs were digested with Rnase R (Epicentre, Madison, WI, USA). A random priming method with the Arraystar Super RNA Labeling Kit was utilized to amplify the enriched circRNAs and to transcribe the circRNAs into fluorescent cRNA. Labeled cRNAs were hybridized onto Mouse circRNA Array V2 (8x15K, Arraystar). Then the slides were washed, and the arrays were then scanned by Agilent Scanner G2505C, and the images were analyzed by an Agilent Feature Extraction software (version 11.0.1.1). Detected circRNAs were regarded as significantly differentially expressed by the value of fold change≥1.5 and P<0.05.

Validation of quantitative real-time reverse transcription-polymerase chain reaction
To validate the data of microarray analysis, quantitative real-time reverse transcription-polymerase chain reaction (qRT-PCR) was performed as previously described [22]. Briefly, total RNAs were transcribed into cDNAs by the SuperScript III Reverse Transcriptase kit (Invitrogen, Carlsbad, CA, USA). The ViiA 7 Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) as well as 2× PCR Master Mix were used to perform the qRT-PCR. The primers used for qRT-PCR were listed in Table 1. The relative expression levels of circRNAs were normalized to GAPDH.

Bioinformatics analysis
According to TargetScan (http://www. targetscan.org) and miRanda (http://www. microrna.org), the miRNA target prediction software (Arraystar) was used to predict circRNA-miRNA interactions. Cytoscape was used to predict a circRNA-miRNA-mRNA network to further visualize the interactions. To forecast the functional annotation of target genes, Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were used for further bioinformatics analysis.

Statistical Analyses
In the present study, statistical differences were assessed by Student t-test, and P<0.05 was considered as statistically significant.

Results
Expression profiling of altered circRNAs in OIR retinas.
CircRNA expressions between OIR retinas and room air control retinas were measured by microarray. The analysis of hierarchical cluster identified different circRNA expression levels in OIR group and control group. A box plot showing the circRNA profiles ( Fig. 1A) indicated similar distributions of all included samples in both groups. The variation of these detected circRNAs between each group was assessed by the scatter plot ( Fig. 1B) and the volcano plot (Fig. 1C). In addition, the hierarchical cluster analysis showed the altered circRNA expression levels as well as the classifications in OIR group and control group. Totally 539 circRNAs were significantly altered in OIR retinas compared with control retinas (fold change≥ 1.5, P<0.05). Of which, 185 circRNAs were significantly up-regulated, and 354 circRNAs were significantly down-regulated in the OIR retinas (Fig.  1D). The top 20 significantly up-and down-regulated circRNAs were listed in Table 2-3. In particular, mmu_circRNA_007438 and mmu_circRNA_012434 were up/down-regulated circRNAs with the most significant alterations in OIR retinas.

Prediction of circRNA-miRNA interactions
CircRNA acted as microRNA sponge associating with related miRNAs, and together they made up circRNA-miRNA axis involving in disease pathogenesis. To determine the function of the validated circRNAs (mmu_circRNA_002573, mmu_circRNA_011180, mmu_circRNA_016108 and mmu_circRNA_22546), their target miRNAs were predicted by informatics analysis based on TargetScan as well as miRanda. The top 5 miRNAs related to each validated circRNA were shown in Fig. 3A, and the predicted interaction sites of mmu_circRNA_002573 and mmu_circRNA_016108 were also displayed in Fig. 3B-C.

Competing endogenous RNA (ceRNA) regulatory networks with GO enrichment and KEGG pathway analyses
To further elucidate the underlying mechanism of altered circRNAs we identified, construction of a circRNA-miRNA-mRNA regulatory network was generated by Cytoscape (Fig. 4). The network was constructed with 236 mRNAs, 4 circRNAs as decoys, and 42 predicted miRNAs in total. The network suggested that circRNAs could indirectly regulate miRNA target genes by competitively binding to miRNA through binding sites.
The GO analysis and KEGG pathway analysis were conducted to explore the biological function of the parental genes. The GO analysis revealed that the term with most genes and the term with the highest enrichment score were both cellular process (GO:0009987) for biological process (Fig.5A). For cellular component, the term "cell part" (GO: 0044464) was with most genes, and was the most significant enriched term as well (Fig.5B). Furthermore, binding (GO:0005488) was related with most genes and the term enriched most was phosphoric ester hydrolase activity (GO:0042578) (Fig.5C) for molecular function. Moreover, the enriched pathways analyzed by KEGG were MAPK signaling pathway, long-term potentiation, renin secretion, spliceosome, pentose phosphate, proteoglycans in cancer, fluid shear stress and atherosclerosis and steroid hormone biosynthesis (Fig.5D). These pathways were related to angiogenesis and endothelial cell motility.

Discussion
Several studies have revealed roles of non-coding RNAs, particularly miRNA and lncRNA [13,23]. Some miRNAs, such as miR-23 and miR-27 served as promotors for neovascularization, while other miRNAs like miR-24 and miR-31 acted as inhibitors [23]. Many miRNAs have been revealed to regulate crucial genes involved in the pathogenesis of ischemic retinopathy. For example, miRNA-126 maintained the integrity of the blood-retina barrier [24] and miR-29a inhibited retinal neovascularization to prevent ROP by down-regulating AGT [25].
However, another spectrum of non-coding RNA, circRNA, was considered to be a by-product of error splicing and was overlooked for a long time. With closed-loop structures, circRNAs are more stable than linear RNAs and play essential roles in the pathogenesis of human diseases [26]. A study reported that 529 circRNAs were abnormally expressed in diabetic retinas in human [15]. Another study identified altered circRNAs in serum samples from patients with or without type 2 diabetes mellitus [14]. Systemic dissection of circRNA profile in OIR model may provide clues to roles of circRNA in pathological neovascularization. A few studies showed an intricate interplay between circRNAs and miRNAs. CircRNAs could competitively bind to miRNAs by miRNA-binding sites and the target genes were regulated as a result [27,28]. We examined circRNA-miRNA interactions and discovered that each selected circRNA, containing at least one miRNA binding sites, was able to interact with several miRNAs (Fig. 3). Our results suggested that circRNAs participated in OIR putatively through targeted miRNA and indirectly regulated gene expression. In addition, there is a subpopulation of circRNAs named exon-intron circRNAs (EIciRNAs), crosslink to RNA polymerase II, that remain in the nucleus and can regulate gene transcription directly in cis [29]. However, the role of EIciRNAs in OIR models largely remains unclear.
The four validated circRNAs by qRT-PCR were selected to draw a whole picture of ceRNA regulatory networks. This circRNA-miRNA-mRNA network may provide clues to the regulatory pathways in OIR for the cascade-amplifying synergistic effects of circRNA-miRNA and miRNA-mRNA. Furthermore, GO enrichment analysis and KEGG pathway analysis were also performed to functionally annotate the predicted target genes. The ceRNA network revealed a novel interaction between the altered circRNAs and 236 mRNAs, and also showed that circRNA can regulate mRNAs through circRNA-miRNA binding. This network served as a shred of credible evidence that circRNA played a role in the pathogenesis of OIR by indirectly targeting certain mRNAs. GO analysis revealed that the most significant enriched term was cellular process in biological process, while cell part was the first highly ranked enriched term in cellular component and phosphoric ester hydrolase activity in molecular function. The KEGG pathway analysis revealed several important pathways related to angiogenesis, including MAPK signaling pathway and renin-angiotensin system (RAS). A number of studies reported that angiogenesis can be restrained via downregulation of MAPK signaling pathway [30][31][32]. Moreover, RAS was an essential growth factor which can stimulate angiogenesis in ischemia [33] and blockage of RAS could inhibit cancer angiogenesis [34]. Altogether, GO and KEGG analysis suggested that mmu_circRNA_002573, mmu_circRNA_011180, mmu_circRNA_016108 and mmu_circRNA_22546, may take part in the process of neovascularization via different pathways.
Although we systemically profiled circRNA expression in retinas of OIR mice, limitations were existed in the study. For example, the circRNA expression profile should be interpreted with caution due to a limited number of samples. Secondly, false-negative data could be presented because the small quantity of circRNAs may below the detection threshold [35]. Moreover, it is possible that the microarray-based screening might miss out some key circRNAs involved in retinal neovascularization, as microarray assays lack the sensitivity of advanced approaches such as next-generation sequencing. Apart from the limitations, there is a promising prospect for circRNAs in diagnosis and therapies since circRNAs are expressed widely and stably. For instance, circ_0005015 promoted retinal angiogenesis by regulating migration, proliferation as well as tube formation of the endothelial cells, and was able to act as potential biomarkers for diagnosis of diabetic retinopathy [15]. A study suggested that exogenous introduction of circRNAs may stimulate the immune system for therapeutic purposes [36]. In addition, circRNA HIPK3 promoted endothelial proliferation and vascular dysfunction by inhibiting miR-30a function [37]. After intravitreal injection of AAV vectors encoding circRNA in mice, multiple retinal cell layers were observed with robust transgene expression, potentially providing novel therapeutic platforms of circRNAs [38].
In conclusion, our study demonstrated that the interactions of circRNAs and miRNAs may indirectly regulate gene expression and altered circRNA may play a role in retinal neovascularization. Further research on downstream molecular functions and their mechanisms may expose potential therapeutic targets for ischemia-induced retinal neovascularization.