Molecular Vision 2017; 23:457-469
Received 26 December 2016 | Accepted 18 July 2017 | Published 20 July 2017
Junde Han, Lingqi Gao, Jing Dong, Jie Bai, Mazhong Zhang, Jijian Zheng
Shanghai Children's Medical Center, Shanghai Jiao Tong University School of Medicine, Shanghai, China
Correspondence to: Jijian Zheng, Department of Anesthesiology, Shanghai Children's Medical Center, Shanghai Jiao Tong University School of Medicine, Shanghai, China 1678 Dongfang Road, Pudong, Shanghai, China, Zip code: 200127; Phone: +86(21) 38626161; FAX: +86(021) 58393915; email: email@example.com
Purpose: Physiologic neuronal apoptosis, which facilitates the developmental maturation of the nervous system, is regulated by neuronal activity and gene expression. Circular RNA (circRNA), a class of non-coding RNA, regulates RNA and protein expression. As the relationship between circRNA and apoptosis is unknown, we explored changes in expression patterns of circRNA during physiologic neuronal apoptosis.
Methods: High-throughput sequencing was used to explore changes in the expression of circRNA in the postnatal developing rat retina. Neuronal apoptosis was determined with immunohistochemistry and terminal deoxynucleotidyl transferase dUTP nick end labeling (TUNEL) in the rat retinal ganglion cell layer.
Results: In total, 2,654, 7,201, and 5,628 circRNA species were detected in the postnatal day (P)3, P7, and P12 rat retina, respectively. Of these circRNA species, 1,371 changed statistically significantly between P3 and P7 and 1,112 changed statistically significantly between P7 and P12. Normal developmental apoptosis, measured with the ratio of apoptotic (caspase-3- or TUNEL-positive) cells to normal cells, showed an increase from P3 to P7 and then a reduction from P7 to P12. In addition, 15 circRNAs whose host genes were associated with apoptosis were differentially expressed during the early development period.
Conclusions: These results associate circRNAs with neuronal apoptosis, providing potential mechanisms and treatment targets for physiologic and drug-induced apoptosis in the developing nervous system.
During the period of synaptic formation and refinement, especially at the peak of synaptogenesis, many neuronal cells undergo apoptosis, or programmed cell death (PCD, also called physiologic apoptosis), to ensure the establishment of accurate neuronal connections and networks [1-4]. In the rodent brain, this process mainly occurs within the 2 weeks after birth, primarily during postnatal days 4 to 14 (P4–P14) [5,6]. The nervous system is also vulnerable to ethanol, general anesthetics, and other substances during this period [7,8]. The underlying mechanism regulating the neuronal apoptotic process remains unclear, but numerous studies demonstrate that neuronal activity and genetic modulation play highly significant roles.
In the developing rat retina, the peak of physiologic apoptosis coincides with the transition from early cholinergic-driven, synchronized spontaneous network activity to glutamatergic-driven activity [9,10]. Blockade of this neural activity by ketamine during this transition period aggravates physiologic apoptosis . During this transition, some relevant receptors, signaling pathways, and apoptosis-related genes also experience dramatic changes [10,12-15]. Previous studies on mRNA or transcription factors to cleave or downregulate mRNA have shown that noncoding microRNA (miRNA) species are involved in neurogenesis, proliferation, axon extension (e.g., mir-9) [16,17], and neuronal apoptosis (e.g., mir-21, -23, -26, -27, and -29) [18-21].
Circular RNA (circRNA), formed by back-splicing, has been reported for decades because of splicing errors [22,23]. More recently, studies have shown that circRNAs are expressed in a variety of eukaryotic organisms, including mammals, and are involved in various physiologic or pathological processes. CircRNAs are widely expressed and show temporal and spatial changes during development [24-27]. This large class of RNA has regulatory abilities, acting as protein or miRNA “sponges” and regulating mRNA transcription or translation [28,29]. Recent studies have shown that circRNAs are remarkably enriched in the nervous system, especially in synapses, during development; circRNAs were also found to regulate synaptic function and neuronal plasticity [30,31]. However, the relationship between circRNAs and physiologic neuronal apoptosis is largely unknown.
Development-related neuronal apoptosis in the central nervous system of rodents mainly occurs within 2 weeks after birth and peaks at about P7 . Thus, we chose to analyze the rat retina at the P3, P7, and P12 time points, to determine any possible associations between physiologic neuronal apoptosis and changes in the expression of circRNA during postnatal development.
Sprague-Dawley rat pups (male and female) aged P3, P7, and P12 days were provided by the Experimental Animal Center at the Shanghai General Hospital in Shanghai, China. Eight rat pups were used for each group of P3, P7, and P12. Rat pups were housed with their dam under a 12 h:12 h light-dark cycle at room temperature 35–37 °C, with food and water available ad libitum until the experiments were conducted. All experimental procedures were reviewed and approved by the Animal Care Committee at the Shanghai General Hospital, Shanghai Jiao Tong University School of Medicine, and were conducted under the guidelines of the Care and Use of Laboratory Animals published by the U.S. National Institutes of Health (National Institutes of Health Publication No. 85–23, revised in 1996) and the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research. Every effort was made to minimize the number of animals used, and their discomfort, during all experimental procedures.
Because of the limited size of the rat retina, we combined the retinas of each group of P3, P7, and P12 rats for high-throughput sequencing of circRNA and immunohistochemistry or terminal deoxynucleotidyl transferase dUTP nick end labeling (TUNEL) assay. The basic experimental protocol was slightly modified from that in previous studies . Briefly, the pups were euthanized instantaneously by decapitation, and then their eyes were rapidly dissected with fine scissors and transferred to an ice-cold (0‒4 °C) bath of artificial cerebrospinal fluid (ACSF: composition in mM: NaCl 119, KCl 2.5, K2HPO4 1.0, CaCl2 2.5, MgCl2 1.3, NaHCO3 26.2, and D-Glucose 11). The bath was continuously bubbled with a 95% O2/5% CO2 gas mixture. About one fifth of the eyeball circumference at the edge of the cornea and sclera was resected using ophthalmological scissors to facilitate the perfusion of ACSF through the retina. After 1 h of recovery in ACSF (37 °C) bubbled with a 95% O2/5% CO2 gas mixture, the retinas were dissected and then either frozen quickly with liquid nitrogen for circRNA high-throughput sequencing or fixed in 4% paraformaldehyde at 4 °C for 24 h for immunohistochemistry and the TUNEL assay.
After the retinas were embedded in paraffin, retinal sections of 4 µm thickness were prepared. Three percent hydrogen peroxide was used to inactivate endogenous peroxidases, and EDTA, adjusted to pH 9.0, was used for heat-induced antigen retrieval for 8–10 min. After rinsing with PBS (1X; 137 mM NaCl, 2.7 mM KCl, 10 mM Na2HPO4, 2 mM KH2PO4, pH 7.4), retinal sections were incubated with rabbit anti-cleaved caspase-3 primary antibody (#9661s, dilution: 1/300, Cell Signaling Technology, Danvers, MA) was performed at 4 °C overnight, followed by incubation with goat anti-rabbit immunoglobulin G (IgG; PV-9001, ZSGB-BIO, Beijing, China) secondary antibody at 37 °C for 1 h. After washing, 3,3′-diaminobenzidine (DAB, ZLI-9017, ZSGB-BIO) was used to visualize immunoreactivity. All sections were then counterstained with hematoxylin to stain the nuclei. Finally, sections were dehydrated and mounted for microscopic examination.
For the TUNEL assay, retinal sections were deparaffinized and rehydrated. After rinsing with PBS, the sections were treated with proteinase K (Roche Applied Science, Indianapolis, IN) and quenched with 3% hydrogen peroxide; then they were incubated in a terminal deoxynucleotidyl transferase (TdT) reaction mix (Roche Applied Science) for 1 h at 37 °C. Sections were then incubated for 5 min with 4',6-diamidino-2-phenylindole (DAPI) at room temperature.
RNA was isolated using the RNeasy system (Qiagen, Dusseldorf, Germany), including on-column DNase I digestion (Qiagen) to eliminate DNA. Then the total RNA was depleted of ribosomal and linear RNA using the RiboMinus kit (Thermo Fisher Scientific, Waltham, MA) and RNase R (Epicenter, Madison, WI), respectively. CircRNAs were validated with quantitative reverse transcription (RT)–PCR using the VAHTSTM Library Quantification Kit (Vazyme Biotech Co., Ltd, Piscataway, NJ). Specifically, samples were heated to 70 °C for denaturing and then cooled to 40 °C on a thermocycler with RNase R buffer to digest linear RNA. After RNase R digestion and circRNA reverse transcription, the cDNA library was produced by PCR (conditions: 98 °C for 30 s, followed by 15 cycles of 98 °C for 10 s, 60 °C for 30 s, and 72 °C for 30 s, then 72 °C for 5 min) following depletion of the second cDNA chain with the USER enzyme. The sequencing libraries were constructed using the VAHTSTM Library Quantification Kit (Vazyme Biotech Co) and were sequenced with an Illumina HiSeq™ 2000 flowcell sequencer (Illumina, San Diego, CA).
The RNA-sequenced (RNA-seq) data were analyzed using FastQC software to verify the quality of the data . The BWA-MEM alignment algorithm of the software BWA (the Burrows-Wheeler Alignment tool) was used to match the RNA-seq data to the genome . To annotate the validated circRNAs, we used CIRI (CircRNA Identifier) software to predict the circRNAs  present in the sample and sequenced them using circBase . The false discovery rate (FDR) was used to correct for multiple hypothesis testing, and a twofold change and FDR≤0.001 was considered statistically significant . Caspase-3-positive or TUNEL-positive cells were counted in a double-blinded manner from five discontinuous views randomly selected from each sample (n = 5 retinas). Data in the figures were analyzed with GraphPad Prism 5 software (GraphPad Software Inc., San Diego, CA) and presented as mean ± standard error of the mean (SEM). A one-way ANOVA (ANOVA) followed by the Bonferroni correction, or the appropriate equivalent non-parametric test, was used for comparisons among groups, and a p-value of < 0.05 was denoted as statistically significant.
We analyzed the raw reads or raw data of the expression of circRNA at P3, P7, and P12, using high-throughput sequencing followed by filtration and quality control (QC; Figure 1A, Table 1) to obtain clean reads as in a previous study . Because exon–exon junction reads (Figure 1B-C) are important features of circRNAs that can map to the junction points of exons, the BWA software was used to match the clean reads to the genome and identify the junction reads. The results showed that the clean reads were expressed abundantly (7–9 × 107/group) and largely mapped to the genome (>75% of all reads) in the early developmental retina (Figure 2A, Table 2). Most of the clean reads mapped on the exon regions could be visualized at P7 (56.41% of all reads) but were present in only small amounts at P3 (15.21% of all reads) and P12 (17.78% of all reads; Figure 2A, Table 2). After the CIRI predictions and false-positive read filtration, validated circRNAs were obtained and annotated according to circBase (Figure 1E, Table 3, Appendix 1, Appendix 2 and Appendix 3).
To identify circRNA changes in the developing rat retina, the total circular junction reads of the different stages were further separated with BWA. CircRNAs with the same sequence were regarded as one circRNA species. The number of circular junction reads at P12 (143,815) was greater than that at P3 (78,758) and P7 (78,385), while the number of circRNA species at P7 (7,021) was greater than that at P3 (2,654) and P12 (5,628; Figure 2B, Table 4).
After the data were integrated, a total of 10,890 circRNAs were obtained; additionally, the overall expression of circRNAs was different among the P3, P7 and P12 time points. For example, there were 4,123 circRNAs expressed solely at P7 (Figure 2C). By comparing the expression of circRNA between the time points, we detected differential expression (greater than a twofold change, FDR≤0.001; Figure 3A,B, Appendix 4 and Appendix 5) of 1,371 circRNAs (out of 8,223 circRNAs) between P3 and P7 and 1,112 circRNAs (out of 10,016 circRNAs) between P7 and P12. We found that 558 circRNAs and 813 circRNAs (out of 1,371 circRNAs) were upregulated and downregulated, respectively, from P3 to P7. The circRNAs upregulated from P3 to P7 included those originating from the genes Akt3, Ccdc6, Gnaz, and Mpped2; the circRNAs downregulated from P3 to P7 originated from the genes Ly6g6f, Melk, Casc5, Smarcc1, and Pank2 (Figure 3C, Appendix 4). We found that 434 or 678 circRNAs (out of 1,112 circRNAs) were upregulated or downregulated, respectively, from P7 to P12. Upregulated circRNAs included those derived from the genes Pde6a, Unc13b, Pde4b, Pi4ka, Mtr, and Dnmbp, while downregulated circRNAs included those derived from the genes Pank2, Smarcc1, Melk, Ly6g6f, Casc5, Dst, and Pak3 (Figure 3D, Appendix 5). After the circRNAs expressed at only P3 or P12 were filtered out, a total of 438 circRNAs were retained for measuring significant changes in the circRNAs across the time points (Appendix 6, Figure 2D). Of these 438 circRNAs, ten were gradually downregulated, and 17 were gradually upregulated from day P3 to P12 (Appendix 7). In addition, 145 circRNAs were present in only small amounts or were undetectable (48 circRNAs) at P7 relative to P3 and P12. In addition, 266 circRNAs of the 438 circRNAs were present at P7 but were absent or diminished at P3 and P12 (172 circRNAs appeared only at P7, Figure 2D).
Most neurons undergoing PCD exhibit the morphological and biochemical hallmarks of apoptosis, such as cell shrinkage, activation of caspases, and DNA fragmentation . Cleaved caspase-3 is produced early in the apoptotic process and is responsible for executing the apoptotic process, making it a strong biomarker of apoptosis . The TUNEL assay detects apoptosis by labeling the fragmented DNA generated during PCD. For enhanced reliability of the results, we used anti-cleaved caspase-3 immunohistochemistry and the TUNEL assay to detect apoptosis in retinas from rats at ages P3, P7, and P12. The cleaved caspase-3-positive-to-negative cell ratio increased from 1.69 ± 0.25% to 3.19 ± 0.42% between P3 and P7 (p = 0.005); it then decreased to 1.69 ± 0.25% between P7 and P12 (p = 0.008; Figure 4A,C). The TUNEL-positive cell ratio increased from 6.75 ± 0.88% to 12.77 ± 1.17% between P3 and P7 (p<0.001) and then decreased to 7.16 ± 0.62% between P7 and P12 (p<0.001; Figure 4B,D). The results confirmed that normal age-dependent physiologic apoptosis of neurons was occurring.
We identified changes in the expression of circRNA that could be linked to normal developmental apoptosis during in the neonatal rat retina. A total of 15 circRNAs from genes associated with apoptosis were obtained after further filtration and screening (Table 5). We found that five circRNAs (from the genes Optn, Thoc1, Rbm5, Ddx19a, and Bnip2) could not be detected, and two circRNAs (from the genes Pias1 and Naa35) were not expressed at P7 but were abundantly expressed at P3 and P12. Three circRNAs (from the genes Rere, Arhgap10, and Birc6) were expressed specifically at P7, and three circRNAs (from the genes Ranbp9, Epha7, and Braf) could be visualized at P7 but were absent or diminished at P3 and P12 (Figure 5).
In this study, we explored physiologic neuronal apoptosis and the expression of circRNAs during early development in the rat retina. Hundreds of circRNAs were abundantly and differentially expressed during the early development of the rat retina. Expression of circRNA followed several patterns: 1) stage-specific expression or suppression, 2) gradual increases or decreases between time points, or 3) stable expression at all stages of early development. We confirmed age-dependent physiologic neuronal apoptosis was occurring in the developing rat retina. In addition, circRNAs originated from genes that were associated with apoptosis (such as Thoc1, Akt3, Rere, Rbm5, Ranbp9, Pias1, Optn, Naa35, Melk, Epha7, Ddx19a, Braf, Bnip2, Birc6, and Arhgap10) changed dramatically in the developing rat retina, suggesting a link between circRNA and physiologic neuronal apoptosis.
Numerous changes in the nervous system occur during the period of peak developmental apoptosis (P4 to P14) [9,10]. For example, the N-Methyl-D-aspartate (NMDA) receptor intracellular signaling pathway switches from the ERK1/2 to the p38 mitogen-activated protein kinase pathway . CircRNA is non-coding RNA mainly composed of exons and is abundant in eukaryotic cells, with highly conserved spatiotemporal specificity [25,26,39]. Several studies have shown that circRNAs are associated with physiologic and pathological processes, such as development of the nervous system [28,40], atherosclerosis [27,41], and cancer [42-44]. The present data indicate that circRNAs may play an important role in the development of the rat retina.
According to these results, apoptosis increased to a peak at or near P7 in the developing rat retinal ganglion cell layer . Although the mechanisms underlying circRNA function are not well-known, studies by Hansen et al. and Memczak et al. showed that circRNA participates in post-transcriptional regulation by acting as a miRNA sponge, interfering with cognate linear isoforms [28,29,42]. In addition, researchers showed that circRNA can bind with Argonaute proteins  or RNA-binding proteins, or can base-pair with other RNAs [45,46], to interfere with the transcription and translation of mRNA. Some circRNAs can even produce proteins .
Through a bioinformatics analysis of circRNAs at different ages of the developing rat retina, the expression of circRNA in 15 circRNAs originating from genes associated with apoptosis was found to change over time. Previous studies have showed large numbers of circRNAs that are endogenous to mammalian cells, and many of these circRNAs are much more stable than linear RNAs [30, 48]; therefore, we believe the observed differences in the expression of transcripts per million (TPM) of the candidate circRNAs were mostly validated rather than simply caused by outliers. As circRNAs can interfere with their cognate linear RNAs, these circRNAs may be linked to physiologic neuronal apoptosis in the developing retina. For example, the Pias1 gene (also known as STAT1) encodes a protein involved in the G1/S transition of the mitotic cell cycle and is a negative regulator of apoptotic processes. The miRNA rno-miR-204-3p is one of the predicted targets of circRNA from Pias1; upregulation of this miRNA induces glioma cell apoptosis . The present study detected low expression of the circRNA from this gene at P7, which may lead to high levels of this miRNA and silencing or cleavage of its cognate mRNA, thus facilitating the apoptotic process. In another example, the Rere gene encodes a member of the atrophin family of arginine-glutamic acid dipeptide repeat-containing proteins; this protein colocalizes with a transcription factor in the nucleus that triggers apoptosis when overexpressed. Notably, overexpression of Rere-derived circRNA may absorb (sponge) rno-miR-128-3p miRNA, which may result in overexpression of the relevant protein (Mapk14, a proapoptotic protein) and subsequent apoptosis . Thus, we believe that the identified circRNAs are potentially important to the processes of physiologic neuronal apoptosis during maturation of the nervous system, acting as microRNA sponges or through other mechanisms mentioned above.
The present study has the following limitations. We did not precisely classify the data or screen other circRNAs that may be related to apoptosis, nor did we verify the candidate circRNAs and their specific effects in physiologic neuronal apoptosis. In addition, we did not analyze the expression patterns of the host genes of the candidate circRNAs, and their effects on microRNA or mRNA, and the expression profiles of the circRNAs in the adult retina. Such studies should be performed in the future. However, we did predict possible functions or pathways through Gene Ontology (GO) enrichment analysis, pathway enrichment analysis, and miRNA-binding target prediction (data not shown), which provide fertile areas for further study.
In conclusion, we confirmed the temporal dependence of developmental physiologic neuronal apoptosis and found 15 circRNAs that had significant differential expression patterns and may be involved in developmental apoptosis. As the nervous system is highly vulnerable to ethanol, general anesthetics, and other substances during this postnatal developmental period, these results may offer potential mechanisms for regulating drug-induced apoptosis in the nervous system.
Appendix 1. Annotation of P3 rat retina circRNAs
Appendix 2. Annotation of P7 rat retina circRNAs
Appendix 3. Annotation of P12 rat retina circRNAs
Appendix 4. Change of circRNAs between P3 and P7 in rat retina.
Appendix 5. Change of circRNAs between P7 and P12 in rat retina.
Appendix 6. Change of circRNAs betweenP3, P7 and P12 in rat retina.
Appendix 7. CircRNAs betweenP3, P7 and P12 in rat retina
This study was supported by the National Natural Science Foundation of China (Beijing, China, Grant No. 81271263 to Jijian Zheng and 81270414 to Mazhong Zhang) and Shanghai Municipal Commission of Health and Family Planning, Key Developing Disciplines (Shanghai, China, Grant No. 2015ZB0106). The authors declare no competing financial interests. Mazhong Zhang (firstname.lastname@example.org) and Jijian Zheng (email@example.com) are co-corresponding authors for this paper.