This article is mentioned in:

Introduction

Retinopathy of prematurity (ROP) is a vasoproliferative disease that primarily occurs in premature babies and in infants with very low birth weights. Retinopathy leads to a series of sequelae, including ametropias, strabismus, amblyopia, complicated cataracts, secondary glaucoma, retinosis, retinal detachment and even blindness (1,2). It has been shown that the administration of high concentrations of oxygen, which is primarily administered to improve the survival of preterm infants, is associated with ROP, even as early as the first wave of ROP (3). Oxygen toxicity caused by dangerously high oxygen saturations (phase 1) followed by hypoxia-mediated vasoproliferation (phase 2) was established and verified in both animal and human studies (4–7). It is well recognized that the hyperoxia-induced suppression of oxygen-regulated angiogenic growth factors, such as erythropoietin (8) and vascular endothelial growth factor (VEGF) (9), leads to both the arrest of retinal vessel growth and the partial loss of existing retinal vessels (phase 1) (10). Subsequently, the increased metabolic requirements and the attenuated vasculature of the developing retina lead to retinal hypoxia, which compensates by stimulating growth factor-induced vasoproliferation (phase 2); finally, this causes fibrous scar formation, retinal detachment, and blindness (1).

At present, intraocular anti-VEGF therapies and several clinical therapeutics, including surgery, cryotherapy and laser photocoagulation, have effectively improved the treatment outcome of ROP (11,12). Furthermore, numerous studies have ascertained the pathophysiological mechanisms underlying ROP, such as structural and functional anomalies, oxidant stress and the roles of growth factors, including VEGF and insulin-like growth factor (12). Nevertheless, knowledge is still lacking about the specific molecular mechanisms, which underlie vasculogenesis and angiogenesis during ROP.

Strategies such as high-throughput techniques and bioinformatics analyses have also been widely used to investigate ROP mechanisms. For example, as early as 2009, Sato et al (13) revealed the gene expression profile of an oxygen-induced retinopathy (OIR) murine model, a valuable platform used to adequately reproduce the vaso-obliteration and neovascularization phases of ROP (5). Thereafter, several bioinformatics predictions and functional studies focused on microRNAs (miRNAs), an important subtype of small noncoding RNAs that are potential critical regulators during retinal vascular angiogenesis in ROP (14–21).

Currently, a novel class of RNA transcripts named circular RNAs (circRNAs) have emerged as potential regulators in a number of different biological processes and various diseases; they primarily serve as miRNA ‘sponges’ that naturally sequester and competitively suppress miRNA activity (22,23). Recently, a few studies have reported the involvement of circRNAs in retinal development and in several other ocular diseases using a high-throughput sequencing analysis (24–26). For example, Han et al (24) identified the circRNA profile during the postnatal development of rat retinas by RNA sequencing and indicated the association between circRNAs and physiological neuronal apoptosis in the developing nervous system. Singh et al (25) analyzed 74 differentially expressed circRNAs from the 12,532 circRNAs identified in the eyes of cystathionine-β-synthase-deficient mice and identified several miRNAs that are specific to each differential circRNA, suggesting a novel regulatory mechanism in hyperhomocysteinemia-induced ocular disorders. Luo et al (26) identified the alanine, aspartate and glutamate metabolic pathway as a target for microRNA (miR)-184 and predicted the expression patterns of circRNAs in the miR-184-transfected human lens epithelial cell line to provide a novel avenue for understanding the role of miR-184 in ocular diseases. Nevertheless, the expression profiles and regulatory networks of circRNAs in ROP remain largely unknown.

In the present study, a high-throughput RNA sequencing strategy and bioinformatics analysis were used to examine the circRNA profiles and the interaction of the circRNA-miRNA-mRNA networks in OIR mouse retinas. The results of the present study predicted a novel regulatory role of circRNA in RNA-RNA coregulation in OIR, which may provide novel insights into the potential mechanisms underlying the development of ROP.

Materials and methods

Animals

C57BL/6J timed-pregnant dams (2.5±0.5 months old; weight 33±0.2 g; n=6) were purchased from Nanjing SIKERUI Biological Technology Co., Ltd. (Nanjing, China). The mice were kept in specific, pathogen-free conditions at the Animal Center of Soochow University (Suzhou, China), housed at 23±2°C and 50±10% humidity with a 12-h light/dark cycle with free access to food and water. The animals were treated in accordance with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research (27). All procedures employed for the animal experiments in the present study were approved by the Animal Care and Use Committee of Soochow University.

Oxygen-induced retinopathy newborn mouse model

A total of 20 neonatal mice (8 female and 12 male; weight 3±0.5 g) were used in the following study and were housed in the same conditions as the pregnant dams. An OIR newborn mouse model of ROP was constructed as described previously (5). Briefly, 7-day-old healthy neonatal mice were randomly divided into the OIR group (n=10) or the control group (n=10). For OIR model construction, 7-day-old mice and their nursing mothers were maintained in 75% oxygen for 5 continuous days. On postnatal day (P) 12, the pups were returned to normoxic air and kept with their mothers for another 5 days. During the entire procedure, the age-matched mouse pups were kept in normoxic air as the control group. On P17, the mouse pups and their mothers were anaesthetized with 10% chloral hydrate (300 mg/kg body weight; Sigma-Aldrich; Merck KGaA, Darmstadt, Germany) via intraperitoneal injection (14) and were sacrificed by cervical dislocation. The retinas of pups were collected for pathology determination and RNA extraction. During the entire experimental process, there were no signs of peritonitis observed in both the newborn pups and the adult mice. The present study was approved by The Ethics Review Board of The Children's Hospital of Soochow University (Suzhou, China).

Retina isolation and vascular staining

On P17, the mouse pups (n=4/group) were anaesthetized with 10% chloral hydrate (300 mg/kg body weight; Sigma-Aldrich; Merck KGaA) via intraperitoneal injection and were sacrificed by cervical dislocation. The intact ocular globe was enucleated and fixed with 4% paraformaldehyde for 1-h at room temperature. After being washed twice in 1X PBS, the retinas were carefully separated from lens and vitreous under a light dissection microscope (Olympus Corporation, Tokyo, Japan; magnification, ×20). The retinas were subsequently subjected to isolectin B4-594 (cat. no. I21413; Alexa Fluor; Molecular probes, Eugene, OR, USA) staining at 4°C overnight, as described previously (6) and were cut into four quadrants and mounted onto a slide. The staining results were determined and imaged under a fluorescence microscope (Olympus Corporation, Tokyo, Japan; magnification, ×40). Retinal whole-mounts were continuously imaged and stitched together using Adobe Photoshop CS5 (Adobe Systems Europe, Ltd., Maidenhead, UK).

RNA library preparation, whole transcriptome RNA and small RNA sequencing

The total RNAs from retinas of mice with OIR and normoxic controls on P17 (n=3/group) were extracted using TRIzol® Reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and the strand-specific RNA-seq libraries were constructed using TruSeq Stranded Total RNA with Ribo-Zero Gold kit (Illumina, Inc., San Diego, CA, USA) according to the manufacturer's instructions as described previously (28). A 1 µg aliquot of each RNA sample was used for library construction. Briefly, after the removal of the ribosomal RNA using the Ribo-Zero rRNA Removal kit (Illumina, Inc.), the RNA was purified by ethanol precipitation and was fragmented. RNA fragmentation was performed using divalent cations at 94°C in 5 µl rRNA Binding Buffer and 5 µl rRNA Removal Mix-Gold. First strand cDNA was synthesized using 8 µl First Strand Synthesis Act D Mix and SuperScript II Reverse Transcriptase at 25°C for 10 min, 42°C for 15 min and 70°C for 15 min. Second strand cDNA was synthesized using 5 µl End Repair Control and 20 µl Second Strand Marking Master Mix. The double strands cDNA was purified by AMPure XP beads (Beckman Coulter, Inc., Brea, CA, USA). After adenylate 3′ ends and adapters ligation, library fragments were purified using the AMPure XP system (Beckman Coulter) for selecting cDNA fragments of preferentially 150–200 bp in length. DNA fragments was enriched through PCR process, and the PCR mixture was incubated in following progress: 1 cycle 98°C for 30 sec; 15 cycles 98°C for 10 sec, 60°C for 30 sec, 72°C for 30 sec; 1 cycle 72°C for 5 min; hold at 10°C. Finally, products were purified (Agencourt AMPure XP, Beckman Coulter) and library quality was assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Inc., Santa Clara, CA, USA). The libraries were sequenced on the Illumina HiSeq X-Ten sequencing platform (OE Biotech, Inc., Shanghai, China, http://www.oebiotech.com/), 150 bp paired-end reads were generated. For small RNA sequencing, the total RNA was ligated sequentially with 3′ and 5′ adaptors, and RT-PCR was performed using the TruSeq™ Small RNA Sample Prep kits (Illumina, Inc.). After isolation using a 6% Tris-borate-EDTA-PAGE gel, the resulting ligation products were sequenced using the same sequencing platform.

Identification and analysis of OIR-responsive differential circRNAs and miRNAs

circRNAs were predicted using CIRI software (v.2.0.3), as described by Memczak et al (22) and Gao et al (29). Briefly, the clean reads were aligned to the reference genome using Bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml); for the unmapped reads, the junctions were picked out using a back-splice algorithm in CIRI software (v.2.0.3). The total circRNAs were blasted against the mouse circRNAs in the circBase database (30) to identify the total novel circRNAs and the total known circRNAs. The expression level of the circRNAs was measured by mapped back-splicing junction reads per million mapped reads. To screen out the differentially expressed circRNAs in OIR retinas, DESeq software (v.1.18.0) (31) was used to normalize the junction read counts of circRNA for each sample, calculate the fold change (FC) and analyze the differential significance between the OIR and control groups. A FC >2 and P<0.05 were selected as significance thresholds.

In addition, the primary functional annotation was performed according to the annotation of circRNA-originated coding genes. Briefly, Gene Ontology (GO) analysis was performed using the DAVID (6.8) online tool (https://david.ncifcrf.gov/) (32) by calculating the number of circRNAs enriched for each GO item. The significance of circRNA enrichment in each GO item was calculated by a hypergeometric distribution test. To explore the pathways possibly regulated by differential circRNAs, a Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to determine the most commonly enriched cellular pathways (33). In addition, miRanda (34) software (v.3.3a) was used to predict the miRNAs that may be sequestered by differential circRNAs that were acting as sponges, and the circRNA-miRNA target interaction network was plotted using the R network package (v.3.6.0, http://www.R-project.org/) (35). Small RNA sequencing was performed and the clean reads were obtained by removing the low-quality reads, the reads with 5′ primer contaminants and poly (A) tails, the reads without the 3′ adapter and insert tag, and the reads shorter than 15 nt and longer than 41 nt. Subsequently, RNA were subjected to the BLAST search against the Rfam v.10.1 (http://www.sanger.ac.uk/software/Rfam) (36) and GenBank (http://www.ncbi.nlm.nih.gov/genbank/). The known miRNAs were identified by aligning against miRBase v.21 databases (http://www.mirbase.org/) (37) and the unannotated small RNAs were analyzed by mirdeep2 (38) to predict novel miRNAs. Differentially expressed miRNAs were identified with a threshold of P<0.05, which was calculated with the DEG algorithm (v.1.18.0).

Prediction of competing endogenous circRNA (ceRNA) and ceRNA networks (ceRNETs) in the OIR model

To examine the potential regulatory roles of circRNAs serving as ceRNA via the circRNA-miRNA-mRNA network in ROP: i) Differential mRNAs and miRNAs from the same OIR samples were separately identified based on whole transcriptome RNA and small RNA sequencing, followed by bioinformatics analysis (Shanghai OE Biotech, Co., Ltd., Shanghai, China); ii) the GO enrichment analysis of the differentially expressed mRNAs was performed; iii) the correlation between miRNA-target (miRNA-mRNA or miRNA-circRNA) was predicted using the miRanda algorithm, with the following parameters: S≥150, ΔG≤-30 kcal/mol and a strict 5′ seed pairing (39); iv) the ceRNA prediction was performed via base sequence and expression value analysis based on the ceRNA_score formula (ceRNA_score = MRE_for_share_miRNA/MRE_for_ circRNA_miRNA); and v). Finally, the ceRNETs were constructed and mapped using Cytoscape software (version 3.2.1, http://cytoscape.org/). Pearson's correlation coefficient and P-values were calculated based on the expression levels of mRNAs, miRNAs and circRNAs, a Pearson's correlation coefficient ≥0.7 and a P-value ≤0.05 were selected as the correlation analysis parameters, and a P-value ≤0.05 was used as the ceRNA analysis parameter.

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

To validate the reliability of high-throughput RNA sequencing and explore the expression trends of circRNA, miRNA and mRNA, RT-qPCR was performed using Q-SYBR green Supermix (Bio-Rad Laboratories, Inc., Hercules, CA, USA) on an Applied Biosystems 7500 FAST Real-Time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). The total RNA from the retinas of mice with OIR and normoxic controls on P17 (n=3/group) were extracted using TRIzol® Reagent (Invitrogen; Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol (40). First-strand cDNA was synthesized from 2 µg RNA per sample using a Transcriptor First Strand cDNA Synthesis kit (Roche Diagnostics, Indianapolis, IN, USA) with the reaction conditions were 10 min at 65°C, 30 min at 55°C and 5 min at 85°C. For RT-qPCR, denaturation was performed at 95°C for 15 min followed by 40 cycles of 95°C for 10 sec, 60°C for 20 sec and 72°C for 30 sec using SuperReal PreMix Plus SYBR Green (Tiangen Biotech, Co., Ltd., Beijing, China). Certain differential circRNAs, miRNAs and mRNAs associated with angiogenesis were determined, and the expression of RNAs were based on the quantification cycle value (Cq) and the relative expression levels were calculated using the 2−ΔΔCq method (41). GAPDH and U6 served as internal standard controls, and all reactions were performed in triplicate. The primers used are shown in Table I.

Table I. Primer sequences for reverse transcription-quantitative polymerase chain reaction.

Table I.

Primer sequences for reverse transcription-quantitative polymerase chain reaction.

Gene symbol	Forward primer, 5′-3′	Reverse primer, 5′-3′	Product length, bp
mmu-miR-203-3p	GTGAAATGTTTAGGACCACTAG	GCTGTCAACGATACGCTACGTAACG	22
mmu_circ_0013414	TACAGACAGTGACAGAGAAGCA	GGCCCATGATTATTCAGAGG	109
mmu_circ_0012698	CGGCCTTATCTGACCTGTAT	GCAACTGCATCTAATGGTCT	106
mmu_circ_0005477	CAAGGTTTCCTCCAAGAGATG	GGTCACTTCGCTTCGTTC	122
mmu_circ_0000668	CCAACCATTGAGCAAATGAAAG	GGAAGAGAGAATGCTGTGGA	108
Srpk2	CATTGCCTTGCTTATTACCCT	CCTTTAGCTTTGTTGACACTAC	81
Rtn4	GGCATCAGGCGCAGATAGAT	CCGTAACCAAGGGAGTGTCC	165
Ptprb	TGGATGTTTGGTTGTACCTACT	GATCTATTACCCTTCACACACG	84
Nrp1	TCAGCTCTGTTTACCAAGCA	ATGGAGGTCACAAACACGA	83
Mapk14	TACTTTGCTCAGTACCACGAC	AGGCTCTTCCACTCATCTAT	101
Kdr	TGGGATGAAATCTTTGGTGGA	ATATCAGGAGCTGGGTAACT	84
Hmox1	CAAGTTCAAACAGCTCTATCGT	GGTCTTAGCCTCTTCTGTCAC	88
Hif1a	TCTGTTCCCATTAGCAGGTGA	AGATTCTGACATGCCACATAG	86
Fn1	ACCACCTCTATCCTCATGT	CAGAACTCTCCTGGAACG	105
Flna	GGTCTCAGTCAAGTTCAACG	ACCCTGACTCCTGAAGACTA	116
Efnb2	GCCAAATCCAGGTTCTAGC	TGCGATCCCTGCGAATAAG	94
Cav1	GTGCTCAGCTATGTAGCATC	CAGTCCAGAGAGGTCTCAC	102
GAPDH	ACCACAGTCCATGCCATCAC	TCCACCACCCTGTTGCTGTA	452
U6	CTCGCTTCGGCAGCACA	AACGCTTCACGAATTTGCGT	94

[i] miR, microRNA; circ, circulating RNA; Srpk2, SRSF protein kinase 2; Rtn4, reticulon 4; Ptprb, protein tyrosine phosphatase receptor type B; Nrp1, neuropilin 1; Mapk14, mitogen-activated protein kinase 14; Kdr, kinase insert domain receptor; Hmox1, heme oxygenase 1; Hif1a, hypoxia-inducible factor 1 subunit alpha; Fn1, fibronectin 1; Flna, filamin alpha; Efnb2, ephrin B2; Cav1, caveolin 1; GAPDH, glyceraldehyde-3-phosphate dehydrogenase; U6, U6 small nuclear RNA.

Statistical analysis

All data are presented as the mean ± standard deviation. A student's t-test was used to analyze the significant difference in the expression of circRNA, miRNA and mRNA between OIR and normoxic control groups and Pearson's rank correlation analysis was used for the correlation analysis in ceRNA networks. The statistical tests were performed using SPSS version 19.0 (IBM, Corp., New York, NY, USA). P<0.05 was considered statistically significant.

Results

Overview of the sequencing results and identification of circRNAs in OIR mouse retinas

OIR mice is a widely accepted model of ROP as the retinas in the model mice are very similar to those of premature infants with regard to the highly immature retinal vessel structure which is hypersensitive to hypoxia damage (42). In the present study, 7-day-old (P7) mouse pups were separately exposed to 75% oxygen or normal air from P7 to P12 to inhibit retinal vessel growth and to cause significant vessel loss. On P12, the mice were subsequently returned to normoxic conditions for 5 days. As expected, both retinal neovascularization and normal vessel regrowth driven by the hypoxic avascular retina manifested on P17, as shown in Fig. 1. For example, obvious pathological retinal angiogenesis, such as uneven and disorganized retinal vessels with obvious dilatation and tortuosity was observed in OIR newborn mice on P17 compared with the control group.

Figure 1. Visualization of OIR in P17 mice. Retinal whole-mounts in the P17 NC and OIR groups are shown by isolectin B4-594 fluorescent dye staining. Inset, the dysregulated neo-vessels in the OIR group were locally magnified. Scale bar, 10 µm in the retinal whole-mounts and scale bar, 1 µm in the magnified inset. NC, negative control; OIR, oxygen-induced retinopathy.

To investigate whether circRNAs was involved in ROP, the extracted RNAs from the mouse retinas were first subjected to high-throughput transcriptome RNA sequencing on P17, and 27,878 candidate circRNAs were identified. Of all the candidates, 18,695 (67.1%) had a length between 200 and 1,000 bp (Fig. 2A), 19,060 (68.4%) had no more than 5 exons (Fig. 2B), and all candidate circRNAs were found to be distributed among all the chromosomes (Fig. 3A). According to the genomic loci of these circRNAs, 25,583 (91.77%) candidates were categorized as sense-overlapping circRNAs, 1,010 (3.62%) were determined as exonic circRNAs, 768 (2.75%) candidates were found to be intergenic, 326 (1.17%) were antisense circRNAs, and 191 (0.69%) were identified as intronic circRNAs (Fig. 3B). In addition, 7,767 (27.9%) circRNAs were identified as known circRNAs by performing a BLAST search against the known mouse circRNAs in the circBase database (Fig. 3C).

Figure 2. Analysis of the genome features of the predicted circRNAs. (A) Length distribution and (B) exon number distribution of the predicted circRNAs were analyzed. circRNA, circulating RNA.

Figure 3. Additional genome features of the predicted circRNAs. (A) Chromosome distribution, (B) subgroup classification and (C) known circRNA identification via circBase of the predicted circRNAs were also analyzed. circRNA, circulating RNA; chr, chromosome.

Functional annotation of differential circRNAs and mRNAs in the OIR model

To identify the possible roles of circRNAs in ROP, the 236 differentially expressed circRNAs, 67 upregulated and 169 downregulated circRNAs, in OIR retinas were screened using the DESeq software (31) (Fig. 4A and Table SI) and were functionally annotated. Using GO analysis, the differential circRNAs were shown to be primarily enriched in the biological processes of the positive regulation of cholesterol esterification, very-low-density lipoprotein particle remodeling and lipoprotein metabolic processes (Fig. 4B and Table II). In addition, using KEGG analysis, several pathways were enriched, such as insulin secretion, Fc gamma R-mediated phagocytosis, lysine degradation, neuroactive ligand-receptor interaction, serotonergic synapse, bacterial invasion of epithelial cells, hepatocellular carcinoma, RNA degradation, non-small-cell lung cancer, Cushing's syndrome, and EGFR tyrosine kinase inhibitor resistance. (Fig. 4C and Table SII). These results provided the potential regulatory cellular pathways by which the differential circRNAs may be involved in the development of ROP.

Figure 4. Functional prediction of differential circRNAs in OIR retinas. (A) Hierarchical cluster of the circRNAs that were differentially expressed in the control and OIR retinas. The log2 value ranges from −2.0 to 2.0 are shown as shades of red or blue according to the key, and P<0.05 was used as the cut-off. Rows represent the individual probe set, and the columns represent the experimental samples. N1-3, sample 1–3 from the control group and O1-3, sample 1–3 from the OIR group. (B) Top 10 upregulated GO terms of all the differential circRNA transcripts listed according to their P-values. (C) Top 20 KEGG predictions for differential circRNAs based on their P-values, and the homologous coding genes from the same genomic loci of the predicted differential circRNAs are shown. The x-axis represents the enrichment score. The size of the bubble indicates the number of enriched circRNAs: The larger the bubble, the greater the number of circRNAs. The variation of colors from purple to red represents the P-values for the enrichment analysis. circRNAs, circulating RNAs; OIR, oxygen-induced retinopathy; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Table II. Differentially expressed circulating RNA transcripts associated with the top 10 upregulated GO terms in oxygen-induced retinopathy retinas.

Table II.

Differentially expressed circulating RNA transcripts associated with the top 10 upregulated GO terms in oxygen-induced retinopathy retinas.

A, Biological process
ID	Term
GO:0010873	Positive regulation of cholesterol esterification
GO:0034372	Very-low-density lipoprotein particle remodeling
GO:0042157	Lipoprotein metabolic process
GO:0033344	Cholesterol efflux
GO:0006637	Acyl-CoA metabolic process
GO:0042271	Susceptibility to natural killer cell mediated cytotoxicity
GO:0042632	Cholesterol homeostasis
GO:0008037	Cell recognition
GO:0043542	Endothelial cell migration
GO:0035329	Hippo signaling
B, Cellular component
ID	Term
GO:0035145	Exon-exon junction complex
GO:0005770	Late endosome
GO:0005769	Early endosome
GO:0030139	Endocytic vesicle
GO:0016592	Mediator complex
GO:0005635	Nuclear envelope
GO:0016605	PML body
GO:0005925	Focal adhesion
GO:0005764	Lysosome
GO:0005604	Basement membrane
C, Molecular function
ID	Term
GO:0060228	Phosphatidylcholine-sterol O-acyltransferase activator activity
GO:0047617	Acyl-CoA hydrolase activity
GO:0017154	Semaphorin receptor activity
GO:0004707	MAP kinase activity
GO:0005102	Signaling receptor binding
GO:0004715	Non-membrane spanning protein tyrosine kinase activity
GO:0004713	Protein tyrosine kinase activity
GO:0003779	Actin binding
GO:0042802	Identical protein binding
GO:0004842	Ubiquitin-protein transferase activity

[i] GO, Gene Ontology.

To fully understand the possible mechanisms and signals involved, 9,756 dysregulated mRNAs were identified using the same cDNA libraries constructed and sequenced above, which included 4,760 upregulated and 4,996 downregulated mRNAs compared to those in the control. The differentially expressed mRNAs with a FC>2 and P<0.05 were selected for GO analysis. As shown in Fig. S1A and Table III, the results indicated that angiogenesis ranked in the top 5 upregulated GO terms that were significantly enriched for biological processes; this coincided with the findings of pathological angiogenesis observed in the OIR mouse retinas (Fig. 1). Therefore, the coregulation network of mRNA-circRNA-miRNA associated with angiogenesis in ROP was examined, and 77 differential mRNAs that were enriched in the biological process of angiogenesis were chosen for the following analysis (Fig. S1B).

Table III. Top 10 upregulated GO terms enriched for biological process of all differential mRNAs in oxygen-induced retinopathy retinas.

Table III.

Top 10 upregulated GO terms enriched for biological process of all differential mRNAs in oxygen-induced retinopathy retinas.

ID	Term
GO:1990744	Primary miRNA methylation
GO:1904714	Regulation of chaperone-mediated autophagy
GO:0046635	Positive regulation of alpha-beta T cell activation
GO:0045870	Positive regulation of single stranded viral RNA replication via double stranded DNA intermediate
GO:0001525	Angiogenesis
GO:0007049	Cell cycle
GO:0016567	Protein ubiquitination
GO:0051301	Cell division
GO:0060271	Cilium assembly
GO:0080009	mRNA methylation

[i] miRNA, microRNA; GO, Gene Ontology.

Prediction and identification of differential miRNAs in the OIR model

In Small RNA sequencing, 14 differentially expressed miRNAs (nine known and five novel miRNAs) were identified. Of these nine known differential miRNAs (Table IV), five miRNAs were upregulated (Fig. 5A), while other four miRNAs including mmu-let-7k, mmu-miR-203-3p, mmu-miR-412-5p and mmu-miR-5128 were downregulated in the OIR samples.

Figure 5. circRNA-miRNA-mRNA coregulatory network was constructed using hypergeometric distribution analysis. (A) Hierarchical cluster analysis of upregulated miRNAs in OIR retinas. log2 values from −2.0 to 2.0 and P<0.05 was used as the cut-off. (B) Heatmap of the differential circRNAs that contain the miRNA-response elements of downregulated miRNAs and those that were co-expressed with angiogenesis-associated differentially expressed mRNAs in OIR retinas are shown. N1-3 represents samples 1–3 from the control group and O1-3 represents samples 1–3 from the OIR group. (C) Prediction of the complementary region of mmu_circ_0000668 for the seeding sequence of mmu-miR-203-3p was performed using miRanda. circRNA, circulating RNA; miRNA/miR, microRNA; OIR, oxygen-induced retinopathy.

Table IV. Differential miRNAs in oxygen-induced retinopathy retinas.

Table IV.

Differential miRNAs in oxygen-induced retinopathy retinas.

A, Downregulated
miRNA_ID	FC value	P-value
mmu-let-7k	0.414673	0.030135a
mmu-miR-203-3p	0.450087	0.03282a
mmu-miR-412-5p	0.516986	0.042373a
mmu-miR-5128	0.122681	1.06×10−5c
B, Upregulated
miRNA_ID	FC value	P-value
mmu-miR-126a-3p	1.754258	0.047684a
mmu-miR-143-3p	2.278029	0.00154b
mmu-miR-146a-5p	2.452372	0.009801b
mmu-miR-210-3p	2.317902	0.02668a
mmu-miR-503-3p	2.733501	0.018975a

a P<0.05

b P<0.01

c P<0.001. miRNA, microRNA; FC, fold change.

Potential regulatory roles of circRNAs serving as ceRNA via circRNA-miRNA-mRNA coregulation in the OIR model

CircRNAs may contain sequences with a high level of homology to linear coding RNAs such as mRNAs. By competing with an mRNA for the same miRNA reaction element (MREs), circRNAs regulate the activity of mRNAs by forming a large ceRNA network. In the present study, 77 upregulated mRNAs that were enriched for the biological process of angiogenesis and 39 circRNAs (Fig. 5B) were separately subjected to the target prediction of four downregulated miRNAs using miRanda (43). For example, mmu_circ_0000668 was predicted to react with mmu-miR-203-3p via two MREs (Fig. 5C). Using a correlation coefficient ≥0.7 and P≤0.05 as a cut-off, correlation analysis was performed for circRNA-mRNA regulatory pair prediction.

A total of 9 circRNAs, including mmu_circ_0005477, mmu_circ_0013414, mmu_circ_0012698, mmu_circ_0000668, mmu_circ_0007799, mmu_circ_0005324, mmu_circ_0008737, mmu_circ_0005768 and mmu_circ_0004158 were predicted as ceRNAs, and the ceRNETs were constructed. As shown in Fig. 6A, 4 ceRNAs (mmu_circ_0005477, mmu_circ_0013414, mmu_circ_0012698 and mmu_circ_0000668) were predicted to regulate the expression of mRNAs, such as SRSF protein kinase 2 (Srpk2), neuropilin 1 (Nrp1), heme oxygenase 1 (Hmox1), hypoxia-inducible factor 1 subunit alpha (Hif1a), filamin alpha (Flna), ephrin B2 (Efnb2), mitogen-activated protein kinase 14 (Mapk14), kinase insert domain receptor (Kdr), fibronectin 1 (Fn1), caveolin 1 (Cav1), reticulon 4 (Rtn4), and protein tyrosine phosphatase receptor type B (Ptprb) by competing with mmu-miR-203-3p.

Figure 6. Validation of the circRNA-miRNA-mRNA coregulatory network. (A) circRNA-miRNA-mRNA network associated with angiogenesis. The blue, green, and red nodes represent the circRNA, miRNA and mRNA, respectively. (B) Validation of the selected circRNAs, miRNAs and mRNAs of the coregulation network was performed using reverse-transcription-quantitative PCR. *P<0.05 vs. respective control. circRNA, circulating RNA; miRNA, microRNA.

In addition, the regulation of fibroblast growth factor receptor (Fgfr)1, Fgfr2, and atypical chemokine receptor 3 by mmu-let-7k may be modulated by mmu_circ_0005324, mmu_circ_0008737, mmu_circ_0005768 and mmu_circ_0004158, which possess the MREs of let-7k. These results suggest that mmu_circ_0005477 may capture mmu-miR-203-3p but also mmu-miR-5128, which target 19 angiogenesis-associated coding genes. Certain RNA transcripts involved in angiogenesis-associated ceRNETs (such as mmu-miR-203-3p-associated) were further verified using RT-qPCR assay as shown in Fig. 6B and were primarily consistent with the sequencing results. In the OIR model, mmu-miR-203-3p was significantly decreased, whereas the expression of mmu_circ_0013414, mmu_circ_0005477 and mmu_circ_0000668 were significantly increased. In addition, the mRNAs expression levels of Ptprb, Nrp1, Kdr, Hmox1, Fn1, Flna and Cav1 were also significantly upregulated (Fig. 6B).

Discussion

At present, ROP is the leading cause of childhood blindness. However, the underlying mechanisms are not fully understood. Dysregulated circRNA expression has been frequently reported in human diseases, including in multiple ocular disorders (22,26,44,45). In the OIR model, a widely utilized mouse model for ROP, neonatal mice were exposed to 75% oxygen between P7 and P12, and the hyperoxia caused vessel regression in the central retina and the cessation of normal radial vessel growth (46). The neonatal mice were returned to normal air between P12 and P17, the period during which vascular regrowth and neovascularization is maximal (6). The OIR and control mouse retinas on P17 were chosen in the present study to determine the expression profiles of circRNAs and mRNAs via high-throughput transcriptome sequencing analysis. A total of 27,878 circRNAs were identified, including 236 differentially expressed circRNAs, as well as 9,756 differentially expressed mRNAs between the OIR and control groups. Of the predicted circRNAs, >70% were <1,000 bp in length, 68.4% had 5 or fewer exons and 27.9% were identified as known circRNAs according to the circBase database. These characteristics matched those of the known characteristics of circRNAs (22), suggesting the reliability of these identified circRNAs. The sequencing results made it possible to determine the regulatory roles of circRNAs in ROP development.

Whole-transcriptome sequencing has been recognized as a more valuable methodology than microarray hybridization as it provides annotation of genes/genomes for the reconstruction of genetic interaction networks that are associated with certain cellular functions, growth/development and biological systems (47). In the present study, 9,756 differentially expressed mRNAs and the 236 circRNAs underwent functional annotation and were found to be enriched in biological processes such as primary miRNA methylation, the regulation of chaperone-mediated autophagy, the positive regulation of α-β T cell activation, the positive regulation of single-stranded viral RNA replication via a double-stranded DNA intermediate and angiogenesis. Angiogenesis ranked in the top 5 upregulated GO terms, in agreement with the well-established development of ROP that is characterized by abnormal angiogenesis (7). During ROP, development of the physiological retinal vasculature is delayed, followed by abnormal retinal neovascularization (7). As a result, pathological angiogenesis leads to retinal fibrosis, detachment and blindness (7,48). These results suggest a possible regulatory mechanism of circRNA, by which, the circRNAs interact with the mRNAs and modulate angiogenesis during ROP. Therefore, a comprehensive understanding of the circRNA-mRNA network may be required to better understand the underlying mechanisms during ROP in addition to the known molecular mechanisms.

Previously, the ceRNA hypothesis has provided a novel angle to interpret RNA-RNA interactions (49). Functional RNA molecules are silenced by miRNA-mediated miRNA-mRNA interactions through sequence-specific binding sites. According to this hypothesis, lncRNAs, circRNAs, pseudogenes, or other mRNAs that contain the same miRNA binding sites are able to competitively bind to miRNA, and, in turn, this regulates the functional activity of the RNAs that are targeted by the now sequestered miRNA molecules (49). As a novel class of noncoding RNAs, circRNAs are usually generated from protein-coding genes by alternative splicing, which can form a covalent ring structure (50), and >80% of circRNAs contain the same sequences as their homologous protein-coding mRNAs or exons, enabling them to serve as ceRNAs by absorbing the miRNAs (49,50).

In the present study, 9 circRNAs were predicted to be ceRNAs that regulate angiogenesis-associated coding genes in ROP by competing with miR-203-3p, miR-5128, and let-7k. Of the three differentially expressed miRNAs, multiple studies have shown that miR-203-3p serves a critical role in the development and progression of a number of different types of tumor (51–55) and may regulate placental angiogenesis through inhibiting the expression of VEGFA and VEGFR2 (56). In ocular diseases, miR-203-3p was reported to serve as a novel regulator of zebrafish retina regeneration by repressing pax6b, which is necessary for progenitor cell expansion (57). Additionally, miR-203-3p is thought to be involved in the retinal differentiation of amniotic epithelial stem cells and umbilical cord blood-derived mesenchymal stem cells, possibly by targeting multiple retina development-relevant genes, such as dickkopf WNT signaling pathway inhibitor 1, cone-rod homeobox, RAR-related orphan receptor beta, neurogenic differentiation 1, neural retina leucine zipper and thyroid hormone receptor beta (58). The results of the present study suggest that the upregulated expression of angiogenesis-associated mRNAs, such as Hmox1 (59), Fn1 (60), Flna (61), Nrp1 (62), Kdr (63), Cav1 (64), Ptprb (65), Srpk2 (66), Hif1a (67), Efnb2 (68) and Mapk14 (69), which were predicted to be targets of miR-203-3p, may be regulated by circ_0005477, circ_0013414, circ_0012698 and circ_0000668 in ROP retinas. The bioinformatics results were in accordance with the previous reports where it was determined that miR-203 targets Hif1a and tenascin C-enriched ECM to regulate glioblastoma aggression (70) and suppresses the expression level of Cav1 mRNA in melanoma cells (71). Liu et al (72) demonstrated that circRNA-ZNF609 silencing ameliorates vascular endothelial dysfunction through a cZNF609/miR-615-5p/MEF2A signaling network, and that circular RNA-HIPK3 mediates retinal vascular dysfunction in diabetes mellitus by acting as an endogenous miR-30a-3p sponge (73). However, the expression profiles and functions of circRNAs in ROP has not been reported, to the best of our knowledge. The present study is the first to identify a novel molecular mechanism that involves the circRNA-miRNA-mRNA interaction in regulating the development of ROP. Hypoxia-induced angiogenesis serves a central role in the pathogenesis of ROP (74), which is also critical for tumor progression, vascularization and metastasis in cancer (75); therefore, ROP and cancer may share common pathogenic mechanisms (76). Considering other time points (e.g. P12 and P21) are also vital for the development of OIR, the present study primarily focused on the differential circRNAs and related regulatory networks involved in angiogenesis during OIR on P17, which may have at least provided novel targets for proliferative retinopathies and possibly for human malignancies.

In conclusion, the present study identified a panel of circRNAs that were differentially expressed in ROP, analyzed the potential biological functions and pathways that may be regulated by these differential circRNAs and predicted the circRNA-miRNA interaction pairs. Furthermore, a circRNA-miRNA-mRNA interaction network was determined. circRNAs may serve a novel regulatory role as a ceRNA, which is involved in the regulation of abnormal angiogenesis during the development of ROP. Nevertheless, detailed functional studies are required to comprehensively validate the results and elucidate the potential molecular mechanisms of competing endogenous circRNAs in ROP.

Supplementary Material

Supporting Data

Acknowledgements

The authors would like to thank Dr Qidong Zu (OE Biotech, Inc., Shanghai, China, http://www.oebiotech.com/) for assistance with the bioinformatics analysis of the RNA sequencing data.

Funding

The present study was funded by The National Natural Science Foundation of China (grant nos. 81500733, 81502500, 81501703, 81302382 and 81602181), The Natural Science Foundation of Jiangsu Province, China (Jiangsu, China; grant nos. BK20150292, BK20150293 and BK20150294) and Suzhou People's Livelihood Science and Technology Project (Suzhou, China; grant no. SS201648).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the author on reasonable request.

Authors' contributions

The present study was conceived by HuZ, KF, YJ and HX. The data were collected by HS, YW, XL, JL and HeZ, analyzed by HS, YW, XL, JL, HeZ and YJ, and interpreted by HuZ, MT, XJ, LZ, YS, YH and HX. The experiments were performed by MT, XJ and LZ. The study was managed and coordinated by HX. Resources were provided by HS, YW, XL, JL, HeZ, YS and YH. The study was supervised by YS, YH and KF. The manuscript was written by HuZ, YJ and HX, and reviewed and edited by YJ and HX. All authors read and approved the final manuscript.

Ethics approval and consent to participate

The present study was approved by The Ethics Review Board of The Children's Hospital of Soochow University (Suzhou, China).

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References