Circular RNA expression profiles in neonatal rats following hypoxic-ischemic brain damage

Introduction

Hypoxic ischemic encephalopathy (HIE) can cause neonatal morbidity and mortality (1). Advances in obstetrical and neonatal care has improved the survival rate of HIE resulting in an increase in the number of infants at risk of permanent neurological deficits (1). Therefore, HIE is still considered a threat to the quality of a child's health and life and further research is required to develop effective neurotherapeutic interventions for HIE.

Circular RNAs (CircRNAs) are a unique type of endogenous non coding RNAs that are formed by back-splicing events via protein coding exons (2,3). They are characterized by conservation, stability, abundance and tissue/developmental stage-specific expression (4,5). Previously, it has been reported that circRNAs have multiple biological functions, which are involved in promoting rolling circle translation, controlling the transcription of parent genes, assisting the formation of alternatively spliced mRNA and acting as microRNA (miRNA/miR) sponges (6-10). Altered levels of specific circRNAs are reported to be associated with numerous human diseases, including cancer, ischemia, stroke, neurodegenerative diseases, heart disease and cartilage degradation (11,12). Recent studies (4,13) performed in adult rats demonstrate that cerebral circRNA expression profiles are significantly altered following stroke and contribute to the stabilization of mRNA expression (14). Few studies reported to date have evaluated circRNA alterations in a hypoxic ischemic (HI) injured brain.

To the best of the authors' knowledge, methods used for the identification of novel disease-associated circRNAs include microarray, reverse transcription-quantitative polymerase chain reaction (RT-qPCR), RNA-seq and capture sequencing (15). To determine whether hypoxic-ischemic brain damage (HIBD) influences the expression of circRNAs to identify novel targets for further studies, alteration of circRNA expression profiles in HIBD rat was investigated in this study using RNA-Seq. The function of circRNA may be associated with the known function of the host genes and therefore gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed to predict potential functions of circRNAs in HIBD. Furthermore, miRNAs regulated by differentially expressed circRNAs were selected for further investigation using RT-qPCR. Using a combination of RNA-seq results, bioinformatics analysis and subsequent RT-qPCR, potential roles of these dysregulated circRNAs in pathophysiological processes in HIBD are proposed.

Materials and methods

The study included 4 procedural modules, including library preparation, high throughput sequencing, bioinformatics analysis and RT-qPCR validation (Fig. 1).

Figure 1 The flowchart of the experimental procedures. The scheme includes four steps. First, sample collection, RNA extraction and library preparation. Then, sequencing with Illumina Hiseq X10 and next, differentially expressed circRNAs and host genes in neonatal brains following hypoxic-ischemic brain injury compared with the sham control brains were identified. Subsequently, bioinformatics analysis, including GO and KEGG pathway enrichment were network performed for analyzing the correlation between circRNAs and miRNAs. Finally, RT-qPCR was utilized to validate the sequencing results. Circ, circular; miRNA, microRNA; RT-qPCR, reverse transcription-quantitative polymerase chain reaction; KEGG, Kyoto encyclopedia of genes and genomes; GO, gene ontology; HIBD, hypoxic-ischemic brain damage; rRNA, ribosomal.

Animals

Sprague-Dawley rats, (age 10 days after birth; 18-20 g, without sex selection) were obtained from the Medical Animal Center of Nanjing Medical University (Nanjing, China). All rats were housed for ≥7 days prior to the study in a temperature (22-25°C) and humidity-controlled (50%) animal facility under a 12 h light/dark cycle. Animals had free access to food and water; however, food was withheld overnight prior to surgery. All animal procedures were approved by the Zhongda Hospital Committee on Animal Research and all the experiments were carried out in accordance with the approved guidelines.

Preparation of HIBD rats

HIBD was produced in rats according to previously described methods (16). Briefly, following anesthetization using ether, rats were subjected to ischemia, with the right common carotid artery permanently double ligated and sliced from the middle. Following a recovery time of 1 h, the animals were exposed to hypoxia (8% O2, 92% N2) for 2.5 h and then returned to their normal living environment. Sham control rats were subjected to isolation and tying ligatures around vessels without occlusion and subsequent ischemia. A total of 24 h following HI insult, the rats were sacrificed using sodium pentobarbital intraperitoneally (40 mg/kg) and their ipsilateral hemispheres, including the cortex and hippocampus, were harvested. The tissues were stored in liquid nitrogen (−196°C) until the following experiment was carried out.

RNA isolation and quality control

Total RNA was extracted from the hippocampus using TRIzol Reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA) according to the manufacturer's protocol. The concentration and quality of RNA was assessed using NanoDrop ND-1000 spectrophotometry (NanoDrop, Wilmington, DE, USA) and RNA integrity was determined using 1% denaturing agarose gel electrophoresis.

Library preparation and RNA-seq

RNA-seq analysis was performed at the Decode Genomics Biotechnology Co., Ltd. (Nanjing, China). The sequencing library was prepared following the removal of ribosomal RNA, according to Illumina® TruSeq® RNA Sample Preparation Guide (Illumina, Inc., San Diego, CA, USA). Ligation of the indexed adaptor was performed following the synthesis of double strand cDNA. Following size selection using Agencourt AMPure XP (Beckman Coulter, Inc., Brea, CA, USA), the libraries were quantified and quality assessed using a Qubit® 2.0 Fluorometer with the Qubit® HSDNA BR assay kit (Invitrogen™; Thermo Fisher Scientific, Inc.) and Agilent bioanalyzer (Agilent Technologies, Inc., Santa Clara, CA, USA), respectively. Sequencing was performed using an Illumina HiSeq X-10 (Illumina, Inc.) with 2×150 bp pair-end technology.

Bioinformatic analyses

Quality control of the sequencing was performed using FastQC version 0.11.2 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Raw data of poor quality (>50% bases having a low-quality threshold <20) and adapter sequences were removed using SOAPnuke, to improve results. The clean data was mapped to rat genomes, provided by Illumina iGenomes (downloaded from http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html) with Tophat2 (version 2.0.7) calling Bowtie2 (version 2.1.0) using default settings. The read numbers mapped to each circRNA were counted. RPM Reads Per Million mapped reads, including TopHat mapping and TopHat-Fusion mapping, were used to calculate the expression level of individual circRNA. Differential expression between circRNAs was assessed using DEGseq algorithm (17). The statistical criteria for selecting aberrant-expressed circRNA was defined using a q-value [the false discovery rate adjusted P-value <0.05 with a fold-change >2.0 or <0.5 (18,19)]. CircRNA annotation was based on the Gencode rat genome (v10) (20). The circRNA donator/acceptor site was intersected to annotate gene regions, including coding RNA, non-coding RNA, intron, antisense and intergenic regions, with the BEDTools suite (v2.16.2) (21). The number of reads aligning to the circRNA-specific head-to-tail junctions was used as a measurement of circRNAs expression.

GO analysis provides a controlled vocabulary to describe genes and gene product attributes in any organism (http://www.geneontology.org). This ontology includes three domains: Biological processes, cellular components and molecular functions. Fisher's exact test was used to detect overlapping data between the differentially expressed list and GO annotation list, except for that expected by chance. Pathway analysis was used to map genes to the KEGG pathways.

To investigate the regulation between circRNAs and mRNAs in the neonatal rats following HIBD, interaction networks were studied among the 6 differentially expressed circRNAs, whose host gene was located on the exon and was identified along with their corresponding target miRNAs. CircRNAs-targeted miRNAs were predicted using miRNA target prediction software based on TargetScan (http://www.targetscan.org/vert_71/) and miRanda (http://www.microrna.org/microrna/home.do). The co-expression network was illustrated using Cytoscape (v2.8.1, https://cytoscape.org/). The associated analyses and drawing were performed by Decode Genomics Biotechnology Co., Ltd. (Nanjing, China).

RT-qPCR

Total RNA was extracted from hippocampus samples using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.) and converted to cDNA using PrimeScript RT-PCR kit (cat. no. RR037A; Takara Bio, Inc., Otsu, Japan), according to the manufacturer's protocol. In brief, the reaction system was incubated at 25°C for 10 min, and then at 42°C for 60 min. PCR was performed in a total reaction volume of 25 µl, containing 12.5 µl SYBR Premix Ex Taq (2X), 2 µl cDNA, 1 µl forward primer (10 µM), 1 µl reverse primer (10 µM), 0.5 µl ROX Reference Dye II (50X) and 8 µl double-distilled water. The amplification procedures were as follows: 10 min at 95°C to initiate denaturation; 40 cycles of 5 sec at 95°C, 30 sec at 63°C and 30 sec at 72°C; and a final extension of 5 min at 72°C. Amplification efficiency was evaluated using standard curve analysis. All samples were normalized using the 2−ΔΔCq method against GAPDH and the experiment was repeated in triplicate (22). The primers of the circRNAs and GAPDH are listed in Table I.

Table I DNA sequence of reverse transcription-quantitative polymerase chain reaction primers used for analysis of circular RNA and mRNA level.

Table I

DNA sequence of reverse transcription-quantitative polymerase chain reaction primers used for analysis of circular RNA and mRNA level.

Name	DNA sequence	Product size (bp)
GAPDH (Rat)	F: 5'-GCTCTCTGCTCCTCCCTGTTCTA-3'	124
	R: 5'-TGGTAACCAGGCGTCCGATA-3'
chr11: 34065887|34070810	F: 5'-CCTCTGCCTTGTTTTGCTGT-3'	146
	R: 5'-GTTGGGGAGGAAAGGGTTTA-3'
chr10: 13931236|13935484	F: 5'-TTTTCCACGATGAACAACC -3'	135
	R: 5'-CAGCAGGTGGCAGCTGTAT-3'
chr1: 200899066|201028171	F: 5'-ATTACCCATGCCCAGATGTT-3'	126
	R: 5'-TGCAGTTTTCAAATGGGTCA-3'
chr13: 31092209|31209056	F: 5'-CTAAACAGGGGCTTCTCAGC-3'	88
	R: 5'-GTCATTGTTCTAAGCATCAGTGG-3'

Statistical analysis

Data analyses were performed using GraphPad Prism 5.0 (GraphPad, Inc., La Jolla, CA, USA), Cluster 3.0 (Human Genome Center, University of Tokyo, Tokyo, Japan) and Cytoscape 2.8.1 (https://cytoscape.org/). Data are presented as the with mean ± standard error of the mean and three repeats were performed for each experiment. Statistical differences were determined using Students t-test and one-way analysis of variance. Multiple comparisons between the groups were performed using the Student-Newman-Kuels method. K-means clustering algorithms were performed to classify the samples according to the selected differently expressed genes. All statistical tests were performed as two-sided tests and P<0.05 was considered to indicate a statistically significant difference.

Results

CircRNAs expression profiles in neonatal rats with HIBD

RNA-Seq was used to investigate circRNAs expression profiles from the hippocampus of HIBD neonatal rats. Hierarchical clustering was performed to demonstrate the circRNAs expression pattern between the HIBD and control groups (Fig. 2A). As presented in the box plot, the distribution of the circRNAs expression profiles across all the samples are similar (Fig. 2B). Statistically significant alterations in circRNAs with between the sham and HIBD groups were identified using a volcano plot (Fig. 2C). Compared with the sham group, 66 circRNAs were demonstrated to be differentially expressed [fold change (FC) >2.0 and P<0.05] following HIBD, including 20 upregulated and 46 downregulated circRNAs. The top 30 altered circRNAs are listed in Table II. The differentially expressed circRNAs following HIBD included 70% exonic and 28% intergenic in origin (Fig. 3A). The distributions of altered circRNAs revealed that the upregulated and downregulated circRNAs were transcribed from all chromosomes, with the exceptions of chr12, chr16 and chr20 (Fig. 3B). The downregulated circRNAs were transcribed from all chromosomes with the exceptions of chr12, chr13 chr16 and chr17 and the upregulated circRNAs were transcribed from most chromosomes (Fig. 3B).

Figure 2 Altered expression profile of circRNAs in the hippocampus following HIBD. (A) Hierarchical clustering analysis of circRNAs that were differentially expressed between the sham and HIBD groups. (B) The distribution of the total expression of circRNAs between the sham and HIBD are nearly the same. (C) Volcano plots demonstrate the differentially expressed circRNAs with statistical significance (FC >2.0 and P<0.05). circRNA, circular RNAs; HIBD, hypoxic-ischemic brain damage.

Figure 3 Features of altered circRNAs. (A) Most (98%) of the circRNAs altered following HIBD are exonic and intergenic. (B) The distribution of altered circRNAs in rat chromosomes. circRNA, circular RNAs; HIBD, hypoxic-ischemic brain damage.

Table II The top 30 differential expressed circRNAs following hypoxic-ischemic brain damage.

Table II

The top 30 differential expressed circRNAs following hypoxic-ischemic brain damage.

circRNA ID	LogFC	P-value	Chromosome	circRNA_type	Gene symbol
chr1: 204959155|204961691	5.7672	0.04349	chr1	Intergenic	-
chr13: 31092209|31209056	5.7672	0.04349	chr13	Exon	Cdh7
chr14: 82168512|82184087	5.7672	0.04349	chr14	Intergenic	-
chr17: 5555491|5571375	5.7672	0.04349	chr17	Exon	Agtpbp1
chrY:1142577|1153312	5.7672	0.04349	chrY	Exon	Kdm6a
chr1: 200899066|201028171	5.9858	0.01412	chr1	Intron	Ate1
chr1: 156552340|156969428	6.1756	0.00824	chr1	Exon	Dlg2
chr1: 235393349|235395737	6.1756	0.00824	chr1	Exon	-
chr1: 195833189|195883848	6.3433	0.00487	chr1	Intergenic	-
chr2: 120517204|120520196	6.3433	0.00487	chr2	Intron	Dnajc19
chr19: 26268232|26303010	6.6295	0.00178	chr19	Intergenic	-
chr2: 177947848|177973555	6.6295	0.00178	chr2	Exon	Rapgef2
chr1:195805819|195840151	6.7538	0.00109	chr1	Intergenic	-
chr1: 195658176|195833285	9.0162	0.00000	chr1	Intergenic	-
chr17:31482295|31493278	9.3864	0.00000	chr17	Intergenic	-
chr11: 83225853|83225993	-9.2371	0.00000	chr11	Intron	Vps8
chr5: 90818736|90827570	-8.9711	0.00000	chr5	Intron	Kdm4c
chr3: 171790813|171790868	-7.1180	0.00018	chr3	Intergenic	-
chr14: 108446849|108456923	-6.9588	0.00027	chr14	Exon	Pus10
chr3: 21894168|21898322	-6.8721	0.00043	chr3	Exon	Strbp
chr10: 59294314|59298784	-6.6812	0.00109	chr10	Exon	Ankfy1
chr11: 34065887|34070810	-6.6812	0.00109	chr11	Exon	Morc3
chr10:74336161|74360926	-6.3371	0.00487	chr10	Exon	Gdpd1
chr8: 79234306|79249597	-6.3371	0.00487	chr8	Exon	Rfx7
chr10: 13931236|13935484	-6.2014	0.00824	chr10	Exon	Pkd1
chr6: 29032002|29074704	-6.2014	0.00824	chr6	Exon	Atad2b
chr3: 112713524|112752649	-6.0515	0.01412	chr3	Exon	Ttbk2
chr5: 128083027|128083901	-6.0515	0.01412	chr5	Intergenic	-
chr5: 148050939|148054556	-6.0515	0.01412	chr5	Intergenic	-
chr9: 53646697|53648288	-6.0515	0.01412	chr9	Exon	Mfsd6

CircRNAs gene symbols GO and pathway analysis

To determine the roles of circRNAs in physiological and pathological processes following HIBD, GO and KEGG pathway analyses were conducted for mRNAs transcribed from the parent genes of altered circRNAs. The results indicated that significant enriched and meaningful GO terms in biological process were ‘positive regulation of guanosine monophosphatase (GTPase) activity’, ‘intracellular signal transduction’ and ‘γ-aminobutyric acid (GABA) signaling pathway’ (Fig. 4A). In terms of molecular function, the majority of the altered circRNA associated mRNAs were identified to be ‘protein binding’, ‘Rab GTPase binding’ and ‘adenosine triphosphate binding’ (Fig. 4B). The cellular component (CC) included enriched CC terms, including ‘postsynaptic density’, ‘synapse’ and ‘cytoplasm’ (Fig. 4C). Gene-enriched KEGG pathway analysis suggested pathways that could be involved during the development of HIBD. The majority of these circRNAs were target genes involved in morphine addiction, nicotine addiction, axon guidance, GABAergic synapse and glutamatergic synapse (Figs. 4D and 5). Glutamatergic synapse and axon guidance have been reported to be involved in neuron injury and regeneration following brain injury in rats (23,24).

Figure 4 GO enrichment and pathway analysis of mRNAs originating from the host genes of altered circRNAs. (A) Major enriched and meaningful GO terms according to biological process. (B) Major enriched and meaningful GO terms according to molecular function. (C) Major enriched and meaningful GO terms according to cellular component. (D) The bar plot demonstrates the significantly enriched pathways. GO, gene ontology; circRNA, circular RNAs.

Figure 5 Kyoto Encyclopedia of Genes and Genomes pathway enrichment for the host gene of differentially expressed circRNAs. The majority of these host genes involved in axon guidance, cGMP-PKG signaling pathway, morphine addiction, Ras signaling pathway, Rap1 signaling pathway and GABAergic synapse. cGMP-PKG, cyclic guanine monophosphate-protein kinase G; GABA, γ-aminobutyric acid.

RT-qPCR of circRNA

To verify the RNA-Seq data, 4 differentially expressed circRNAs were selected for qPCR and compared with the sham group, including circRNAs chr13: 31092209|31209056 and chr1: 200899066|201028171, which were upregulated, and chr11: 34065887|34070810 and chr10: 13931236|13935484, which were downregulated significantly following HIBD (P<0.01; Fig. 6). The expression patterns of these circRNAs were consistent with the RNA-Seq data.

Figure 6 QPCR validation results of four selected circRNAs. A total of four circRNAs (chr1: 200899066|201028171, chr10: 13931236|13935484, chr11: 34065887|34070810 and chr13: 31092209|31209056) were selected for validation by qPCR in eight samples. Data are presented as the mean ± standard error of the mean (n=4 per group). *P<0.05 and **P<0.01 vs. the sham group (Student's t test). QPCR, quantitative polymerase chain reaction; circRNA, circular RNA; HIBD, hypoxia induced brain damage.

Annotation for circRNA/microRNA interaction

circRNAs serve a crucial role in sponging and regulating the expression of miRNA, their interaction with HIBD associated miRNAs suggests that circRNAs could serve an important role in pathophysiological processes of HIBD. Therefore, potential circRNA/microRNA interactions were investigated based on TargetScan and miRanda. The circRNA/miRNA interaction of the differentially expressed circRNAs is presented in Fig. 7. These results are in agreement with previously reported differential expression of miRNA in neonatal rats (25-29).

Figure 7 Annotation for circRNA/miRNA interaction. The orange nodes represent the differentially expressed miRNAs and the green nodes represent parts of the circRNAs identified to be differentially expressed in this study. Circ, circular; miRNA, microRNA.

Discussion

To the best of the authors' knowledge, this is the first study to profile circRNA expression in a rat hippocampus following HIBD. Initially, 66 circRNAs were identified to be differentially expressed in the early stages following HIBD. Furthermore, significantly altered GO terms were identified, mostly correlated pathways and predicted circRNA/miRNA interactions using bioinformatics analysis. A total of four circRNAs were verified using RT-qPCR, which was consistent with the results from RNA-seq data. The potential roles of the altered circRNAs in the pathophysiological processes following HIBD suggest they could be potential biomarkers or novel therapeutic targets in the treatment of HIBD.

The reported evidence indicates that >98% transcriptional output is not used in protein translation, but forms various classes of non-coding (nc)RNAs, including miRNAs, long non-coding RNAs and circRNAs (30). It is currently under debate whether ncRNAs have multiple functions, including controlling translation, transcription, epigenetics and RNA/protein scaffolding. In general, ncRNAs are considered to have novel roles in genomic regulation and dysregulation of ncRNAs is associated with the onset of disease and secondary neural damage following central nervous system (CNS) injury. Notably, several pathophysiological mechanisms, including inflammation, blood-brain damage, apoptosis, autophagy, oxidative stress and endoplasmic reticulum stress are closely associated with the expression of ncRNAs, it is suggested that the brain damage may be corrected by determining the level of dysregulated ncRNAs (31-34). Therefore, reversing the altered ncRNAs could be a potential therapeutic strategy for HIBD.

miR-126 was significantly upregulated following HIBD, which was in agreement with our previously report; miR-126 expression levels were associated with retinopathy of premature neonatal rats induced by acidosis, in which mRNA and protein expression levels of the vascular endothelial growth factor were notably upregulated (35,36). Similar results were identified for miR-9a, which was reported to be significantly increased in the plasma of hyperoxia-induced neonatal rats and involved in ‘in utero embryonic development’, ‘neurotrophin signaling pathway’, ‘B cell receptor signaling pathway’ and ‘mitogen activated protein kinase signaling pathway’. On the contrary, miR-21a and miR-25a were downregulated in the injured cortex following hypoxic-ischemic brain damage and contribute to neuronal death with increased expression levels of the proapoptotic B-cell lymphoma 2 family members-Noxa and Bax. In contrast, increasing the expression levels of these miRNAs significantly mitigates neural damage (37,38).

Overall, the results from this study suggest potential in the treatment of HIBD using targeted miRNAs and lncRNAs. However, the role of circRNAs, a new star of ncRNAs, in the pathophysiological process of HIBD is yet to be reported. The expression profiles of circRNA in this study were altered following HIBD. As previous studies report, circRNAs can regulate the transcription of their host genes (39,40) and dysregulation could influence numerous molecular events essential to the process of brain damage following HIBD. Parent genes that generate HIBD-altered circRNAs primarily participated in the positive regulation of guanosine triphosphatase activity, intracellular signal transduction and GABA signaling pathway, which are indispensable for neural regeneration. Synapse, neuron part and neuron projection are essential for neuron-neuron communication, producing movement, feeling and memory. Altered circRNA-associated mRNAs were primarily located in the synapse, and the neuron to neuron projection, which may be correlated to post-HIBD pathophysiology. However, further investigation is required to determine the regulation of circRNAs in HIBD.

CircRNAs are reported to act as miRNA sponges and suppress miRNA activity, resulting in upregulation of miRNA targets (5,8,41). miRNAs serves key roles in normal CNS development and function (42,43) and certain circRNAs could be involved in HIBD through circRNA/miRNA interactions. The dysregulated circRNAs following HIBD was demonstrated to include several miRNA binding sites using miRNA target prediction software. The correlation of miRNAs with HIBD suggest that circRNAs could serve important roles in HIBD. For example, chr1: 200899066|201028171 was predicted to bind miR-126a, miR-9a and miR-26a and miR-9a was identified to be associated with hypoxia-induced neuronal apoptosis (26,44-46).

In conclusion, the results in this study indicated that circRNAs were significantly altered in the hippocampus following HIBD compared with the sham control. Using bioinformatics analysis and circRNA/miRNA interaction prediction, circRNAs could be involved in brain damage and also neural regeneration following HIBD. Therefore, these findings suggest a potential treatment of HIBD through the modulation of circRNAs. Deciphering the biological functions of the circRNAs identified in this study requires further investigation.

Funding

The present study was financially supported by the Project of Clinical Advanced Techniques, the Primary Research and Development Plan of Jiangsu Province (grant no. BE2017719), the Pediatric Medical Innovation Team of Jiangsu Province (grant no. CXTDA2017022), the Project of National Youth Fund (grant no. 81601355) and the Project of Postdoctoral Fund of Jiangsu Province (grant no. 1701162C).

Availability of data and materials

The datasets used in this study are available from the corresponding author on reasonable request.

Authors' contributions

LJ, RZ and ZX contributed to the concept and design of the present study, prepared the manuscript and conducted the experiments and data analysis. HL and RZ performed animal treatments. RZ contributed to the acquisition of data. LJ and ZF conducted bioinformatic analysis. ZX also provided guidances. All the authors read and approved the final version of the manuscript.

Ethics approval and consent to participate

All animal procedures were approved by the Zhongda Hospital Committee on Animal Research and all the experiments were carried out in accordance with the approved guidelines.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that there are no known competing interests regarding the publication of this study.

Acknowledgments

Not applicable.

References