Circular RNA Expression Alteration and Bioinformatics Analysis in Rats After Traumatic Spinal Cord Injury

Spinal cord injury (SCI) is mostly caused by trauma. As primary mechanical injury is unavoidable in SCI, a focus on the pathophysiology and underlying molecular mechanisms of SCI-induced secondary injury is necessary to develop promising treatments for SCI patients. Circular RNAs (circRNAs) are associated with various diseases. Nevertheless, studies to date have not yet determined the functional roles of circRNAs in traumatic SCI. We examined circRNA expression profiles in the contused spinal cords of rats using microarray and quantitative reverse transcription-PCR (qRT-PCR) then predict their potential roles in post-SCI pathophysiology with bioinformatics. We found a total of 1676 differentially expressed circRNAs (fold change ≥ 2.0; P < 0.05) in spinal cord 3 days after contusion using circRNA microarray; 1261 circRNAs were significantly downregulated, whereas the remaining 415 were significantly upregulated. Then, five selected circRNAs, namely, rno_circRNA_005342, rno_circRNA_015513, rno_circRNA_002948, rno_circRNA_006096, and rno_circRNA_013017 were all significantly downregulated in the SCI group after verification by qRT-PCR, demonstrating a similar expression pattern in both microarray and PCR data. The next section of the study was concerned with the prediction of circRNA/miRNA/mRNA interactions using bioinformatics analysis. In the final part of the study, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes analyses indicated carbohydrate metabolic process was one of the most significant enrichments and meaningful terms after GO analysis, and the top two signaling pathways affected by the circRNAs-miRNAs axes were the AMP-activated protein kinase signaling pathway and the peroxisome related pathway. In summary, this study showed an altered circRNA expression pattern that may be involved in physiological and pathological processes in rats after traumatic SCI, providing deep insights into numerous possibilities for SCI treatment targets by regulating circRNAs.


INTRODUCTION
Spinal cord injury is a growing public health concern worldwide that is accompanied by permanent neurological impairment and attendant social and economic losses (Yang et al., 2014;Manohar et al., 2017). SCI is mostly caused by trauma , which involves primary mechanical injury due to rapid direct compression and contusion on the cord and secondary injury, such as hemorrhage, edema, ischemia, cell death and oxidative stress, during the following hours (Zhu et al., 2017b). It is well-known that primary mechanical injury is hardly avoidable. As a result, concentrating on the pathophysiology and underlying molecular mechanisms of SCI-induced secondary injury is necessary to develop promising diagnosis and treatment approaches for SCI patients . However, no neuroprotective and regenerative therapies are currently available that can directly produce beneficial effects. Although past decades have seen the rapid development of stem cell transplantation (Nagoshi and Okano, 2017), surgical decompression (Yang et al., 2013;Hu et al., 2015) and high-dose methylprednisolone (Hextrum and Bennett, 2018) in the field of SCI, a multitude of controversies on the efficacy of these approaches remains.
Recent advances in ncRNAs, including lncRNA and miRNA, may be at the point of breaking this impasse. Increasing studies in this field have examined an outstanding outcome that neurologic damage can be alleviated by normalizing the expression levels of certain lncRNAs and miRNAs in rats after SCI (Hwang et al., 2016;Lv, 2017;Zhu et al., 2017a;Zhou et al., 2018a). Circular RNAs (circRNAs) are a special type of endogenous non-coding RNAs formed by back-splicing events via protein-coding exons (Qu et al., 2015). Notably, they have been attracting a wealth of interest because growing research has found that the altered expression of specific circRNAs was closely related to human diseases, such as cancer (Zhou et al., 2018b), neurological diseases (Bai et al., 2018), and cardiac diseases (Salgado-Somoza et al., 2017). In recent years, a selection of studies stated that circRNA expression profiles were significantly altered after traumatic brain injury in mice and rats (Xie et al., 2018;Zhao et al., 2018). Nevertheless, there are no data concerning circRNAs in traumatic SCI so far, and their molecular and intermolecular interactions and key signaling pathways remain to be elucidated. Although the physiological functions of circRNAs are largely unknown, they are expected to regulate the transcription of parent genes, promote rolling circle translation, and help form alternatively spliced mRNAs and sponge miRNAs, namely, ceRNAs (Memczak et al., 2013;Liu et al., 2017b).
We used high-throughput microarray analysis to screen circRNAs expression patterns in the spinal cord of adult rats after traumatic SCI to determine whether the expression levels of circRNAs were altered and lay a foundation for future work. Subsequently, we selected five differentially expressed circRNAs for low-throughput validation by qRT-PCR. Then, we used bioinformatics tools to predict their putative biological functions with miRNA and mRNA by establishing circRNA/miRNA/mRNA networks and to identify feasible functions of all mRNAs regulated by these specific circRNAs in traumatic SCI.
The purpose of this investigation was to explore the circRNA expression profiles in rats after traumatic SCI, to determine the potential roles of these differentially expressed circRNAs in post-SCI pathophysiology processes.

Animals and Experimental Groups
This study was carried out in accordance with the principles of the Basel Declaration and recommendations of the National Institute of Health Guide for the Care and Use of Laboratory Animals (NIH Publications No. 8023, revised 1978). The protocol was approved by the Institutional Animal Care and Use Committee of Capital Medical University. Female Sprague-Dawley rats (age: 10-12 weeks; weight: 250-300 g) were purchased from the Animal Care Center of Academy of Military Medical Sciences (Beijing, China). Rats were housed in a temperature-controlled (20-28 • C) and light-controlled (12-h light/dark cycle) room. They were habituated to the housing conditions for at least 7 days before SCI. Additionally, the animals had free access to standard rat chow and tap water; however, food was withheld overnight before surgery.
Twelve rats were randomly assigned to two groups using a computer-generated randomization schedule: rats in the sham control group (n = 6) were treated with laminectomy alone without contusion, and rats in the SCI group (n = 6) were subjected to laminectomy plus contusion. The operators performing the surgeries were blinded to the experimental groups. Before and after surgery, our two experienced researchers performed a behavioral test of rats in the two groups using the BBB (Basso, Beattie, and Bresnahan) score, which is frequently used to observe the locomotor function of rats (Scheff et al., 2002). The BBB score ranges from 0 (no hindlimb movement) to 21 (normal movement), and the final score was given by consensus.

Contusion SCI Model and Tissue Collection
In this study, the contusion injury method was adopted to trigger moderate injuries in rat models. Briefly, all rats were anesthetized by intraperitoneal injection of 0.4 mg/g body weight 10% chloral hydrate (Kermel, Tianjing, China), and their backs were shaved and sterilized. After suspending the rats in a stereotaxic frame, a 4-cm-long longitudinal midline incision was made to expose the T9-T11 spinal column. Following stripping of the paraspinal muscles, laminectomy was performed at the T10 level to expose the spinal dura mater without tearing it in a sterile condition. Subsequently, rats were clamped by their spinous processes at T9 and T11 with sterilizing forceps, which was followed by spinal cord contusions induced by an Infinite Horizon Impactor (IH-0400 Impactor, Precision Systems and Instrumentation, LLC, United States), leading to a moderate injury at the T10 level. For parameter setting, a standard rat tip impactor size (2.5 mm in diameter), programmable dwell time (1 s) and programmable force levels (225 kDynes) were applied to induce moderate intensity injuries with the NYU weight-drop device Frontiers in Molecular Neuroscience | www.frontiersin.org (2.5 mm in diameter, 10 g rod, height of 12.5 mm) according to the conversion equation (Basso et al., 1996;Khuyagbaatar et al., 2015). Errors greater than 3% in terms of force levels and fracture of spinal dura mater were not accepted.
After inducing contusions, the operative region was washed gently with warm 0.9% normal saline twice (37 • C, 2 mL) to avoid local infection. The cord surface showed signs of subarachnoid hematoma and an intense dark brown/purple color. The rats with swinging tails that quickly retracted their lower limbs immediately after SCI were regarded as eligible as described previously (Hu et al., 2015). Finally, the wound was then sutured in layers. A warm environment was established to maintain body temperature during surgery. Notably, 4 mL of Ringer lactate solution was administered intraperitoneally to supplement electrolytes and body fluids. Rats were then housed in individual cages. Animals were fed with free access to food and water. Penicillin (40000 U, intramuscular injection) was administered daily for 3 days to prevent systemic infection. The bladders were emptied manually every 8 h until the rats were killed.
At 3 days post-surgery, rats in all groups were euthanized with an overdose of 10% chloral hydrate (10 mL/kg), and a 1-cm long segment of spinal cord, including the injury epicenter, was quickly dissected and collected without transcardial perfusion in advance and was then fresh frozen in liquid nitrogen to prevent RNA degradation similarly to a previously reported method (Pang et al., 2016).

RNA Extraction and Quality Control
Briefly, total RNA was isolated with TRIzol reagent (Invitrogen, Carlsbad, CA, United States) based on the manufacturer's instructions. After identifying the purity and concentration of total RNA from each sample with a NanoDrop ND-1000 (NanoDrop, Wilmington, DE, United States), the integrity of RNA and gDNA contamination was tested using denaturing agarose gel electrophoresis.

CircRNA Microarray
Cord tissue samples from rats in the SCI (n = 3) and sham groups (n = 3) were used in the microarray experiment. As shown in Supplementary Figure 1A, the sample preparations, sample labeling, and microarray hybridizations were performed based on the manufacturer's standard protocols (Arraystar). First, total RNAs were digested with Rnase R (Epicentre, Inc.) to remove linear RNAs and enrich circRNAs. Second, the enriched circRNAs were amplified and transcribed into fluorescent cRNA using a random priming method (Arraystar Super RNA Labeling Kit; Arraystar). Third, the labeled cRNAs were purified by a RNeasy Mini Kit according to the manufacturer's instructions (Qiagen). Then, the concentration and specific activity of the labeled cRNAs (pmol Cy3/µg cRNA) were accessed by a NanoDrop ND-1000. Afterward, 1 µg of each labeled cRNA was fragmented by adding 5 µL of 10× blocking agent and 1 µL of 25× fragmentation buffer, and then, the mixture was heated at 60 • C for 30 min; finally, 25 µL of 2× hybridization buffer was added to dilute the labeled cRNA. Approximately 50 µL of hybridization solution was dispensed into the gasket slide and assembled onto the circRNA expression microarray slide. Lastly, the slides were incubated for 17 h at 65 • C in an Agilent Hybridization Oven. The hybridized arrays were washed, fixed and scanned with an Agilent Scanner G2505C. The microarray experiments in our study were performed by KangChen Bio-tech, Shanghai, China. This standard procedure was also used in other publications (Xie et al., 2018).

CircRNA Data Analysis
Acquired array images were analyzed using Agilent Feature Extraction software (version 11.0.1.1). As shown in Supplementary Figure 1B, quantile normalization and subsequent data processing were performed with the R software limma package. After quantile normalization of the raw data, low intensity filtering was performed. When comparing two groups for profile differences (SCI versus sham control), the "fold change" (i.e., the ratio of the group averages) between the groups of each circRNA was computed. The statistical significance of the difference was conveniently estimated by a t-test. False discovery rate was calculated from Benjamini Hochberg FDR to correct the P.
Differentially expressed circRNAs with statistical significance (fold changes ≥ 2.0, P < 0.05) between the groups were identified through volcano plot filtering or fold change filtering. Hierarchical clustering was used to show distinguishable circRNAs expression patterns among tissue samples using MeV (Multiple Experiment Viewer 1 ). Additionally, the distribution of differentially expressed circRNAs in rat chromosomes and the bar diagram of circRNA categories were performed by GraphPad prism 6.
Next, we first predicted the potential sponging miRNAs for each differentially expressed circRNA with Arraystar's homemade miRNA target prediction software according to TargetScan (Enright et al., 2003) and miRanda (Pasquinelli, 2012) in order to investigate the functions of dysregulated circRNAs. We ranked miRNA candidates based on the mirSVR and listed the five highest-ranking miRNA candidates for each circRNA.

qRT-PCR Assay
The expression of five circRNAs randomly selected from microarray analysis was validated using a qRT-PCR assay. Cord tissue samples from rats in the SCI (n = 6) and sham groups (n = 6) were used in the PCR experiment. As previously reported (Xie et al., 2018), total RNA was isolated using TRIzol reagent according to the manufacturer's standard protocols (Invitrogen). RNA quantification and quality were accessed by a NanoDrop ND-1000, and RNA integrity was identified by electrophoresis on a denaturing agarose gel. Next, we synthesized cDNA in line with the manufacturer's instructions. qRT-PCR was performed in a ViiA 7 Real-time PCR System (Applied Biosystems) with PCR master mix (2×, Arraystar). The parameter settings were 95 • C denaturation (10 min), 95 • C (10 s), and 60 • C (60 s), which was repeated for 40 cycles. After the amplification reaction was finished, the procedure was performed as follows: 95 • C (10 s), 60 • C (60 s), and 95 • C (15 s). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) acted as an internal control to normalize the data. The relative expression levels of circRNAs were calculated using the relative standard curve method (Larionov et al., 2005).

Prediction for CircRNA/miRNA/mRNA Interactions
The ceRNA hypothesis is that RNA transcripts can crosstalk by competing for common miRNAs, where MREs are the foundation of this interaction (Salmena et al., 2011). Any RNA transcript with MREs might act as a ceRNA, and ceRNAs include pseudogene transcripts, circRNAs and mRNAs; these transcripts can compete for the same MER to mutually regulate it. We constructed the ceRNA network by merging common targeted miRNAs. The circRNA/miRNA/mRNA interactions were predicted in combination with in-house miRNA target prediction software according to TargetScan and miRanda software (Friedman et al., 2009;Huang et al., 2016) to identify the potential targets of miRNAs. Three conditions must exist for the ceRNA network to occur (Salmena et al., 2011). First, the relative concentration of ceRNAs and their miRNAs is clearly important; second, the effectiveness of a ceRNA depends on the number of miRNAs that it can "sponge"; third, not all of the MREs on ceRNAs are equal. Therefore, we accept only the ceRNA-pair relations that passed filtering measures (P < 0.05). A graph of the circRNA/miRNA/mRNA network was made using Cytoscape software (version 3.5.1) to visualize these relationships.
Next, we used Gene Ontology (GO 2 ) to reveal the biological process, cellular component and molecular function of the target mRNAs. Significant pathways were identified using pathways in the Kyoto Encyclopedia of Genes and Genomes database (KEGG 3 ). A P < 0.05 indicates the significance of GO and KEGG pathway terms. The FDR was calculated to correct the P.

Statistical Analysis
Statistical analysis was performed using SPSS software (version 21.0, Chicago, IL, United States). The results were shown as the mean ± SEM. Student's t-test was used to compare significant differences between the two groups. A P < 0.05 was considered to be statistically significant.

Database and Accession Numbers
CircRNA microarray data in our study were deposited at the NCBI Gene Expression Omnibus (GEO) under the Accession No. GSE114426.

Expression Pattern of CircRNAs in Spinal
Cord Tissues After SCI anesthesia recovery, which indicates a loss of locomotive function due to the contusion (data not shown).
The first set of questions aimed to identify circRNAs expression levels in the two groups, so a high-throughput microarray assay was used. It can be seen from the box plot ( Supplementary Figure 2A) that the normalized intensity values in all samples were not significantly different, suggesting a similar distribution of circRNAs expression profiles in each group. Supplementary Figure 2B showed RNA integrity and the gDNA contamination test. After microarray scanning and normalization, the dysregulated expression of all 13279 circRNAs between groups was detected. A total of 1676 circRNAs were differentially expressed (fold change ≥ 2; P < 0.05) between the groups, 1261 of which were significantly downregulated and 415 of which were significantly upregulated in the SCI group. Then, Table 1 provides experimental data on the top 10 upregulated and downregulated circRNAs ranked by fold changes after microarray analysis. Two graphs of hierarchical clustering of circRNAs expression were performed (shown in Figure 1) to visualize these data. Next, volcano plots were constructed using fold change values and P values to visualize differential expression between the two different conditions (Figure 2A), and scatter plots presented the circRNA expression variation (or reproducibility) between the SCI and sham control groups ( Figure 2B). Additionally, the distribution of differentially expressed circRNAs in chromosomes is presented in Figure 2C. What stands out in Figure 2D is that the bar diagram of the circRNA category based on gene sources revealed 325 exonic, 6 antisense, 7 intronic, 9 intergenic, and 68 sense overlapping circRNAs that were upregulated. In contrast, 1012 exonic, 4 antisense, 5 intronic, 33 intergenic, and 207 sense overlapping circRNAs were downregulated.
Lastly, we predicted the five most likely potential target miRNAs for each differentially expressed circRNA. The five highest-ranking miRNA candidates as binding targets of each circRNA are listed in Table 2 (for the top 10 upregulated and downregulated circRNAs).

Verification of qRT-PCR
Overall, 1261 differentially expressed circRNAs were listed after we filtered the high-throughput microarray assay data; therefore, the accuracy of the data needed verification. Next, five circRNAs with relatively high fold changes and similar tissue sample distributions within each group were randomly selected for low-throughput verification with qRT-PCR. Compared with the sham control, rno_circRNA_005342, rno_circRNA_015513, rno_ circRNA_002948, rno_circRNA_006096, and rno_circRNA_ 013017 in the SCI group were all significantly downregulated after SCI, suggesting a similar expression pattern in these candidate circRNAs in both the microarray and PCR data, as shown in Figures 3A-E. Table 3 provides the fold change and P-value for each candidate circRNA between groups in terms of microarray and PCR analysis to make data interpretation more straightforward. As presented in Table 4, the sequences of the primers used for qRT-PCR analysis of circRNA are listed. We also provided amplification plots and melt curve plots of these candidate circRNAs and GAPDH in Supplementary Figures 3, 4.    FIGURE 3 | Validation of five selected circRNAs using qRT-PCR. Compared with the sham control, rno_circRNA_005342, rno_circRNA_015513, rno_circRNA_002948, rno_circRNA_006096, and rno_circRNA_013017 in the SCI group were all significantly downregulated after SCI after validation by PCR assay in 12 samples (A-E). The data were normalized using the mean ± SEM (n = 6 per group). * P < 0.05, * * P < 0.01, * * * * P < 0.0001.

Bioinformatic Predictions of CircRNA/miRNA/mRNA Interactions
The five candidate circRNAs validated by qRT-PCR were next selected to construct circRNA/mRNA of miRNA interactions using bioinformatics tools. Notably, an entire network of these interactions was depicted (shown in Figure 4). The results showed that a total of 5 circRNAs, 60 miRNAs and 253 mRNAs were included, presenting a large interaction network. CircRNAs

GO and Pathway Analysis of Putative Target mRNAs
Based on these interactions, five candidate circRNAs may play an important role in molecular mechanisms by regulating target miRNAs/mRNAs. As a result, we performed GO and pathway analysis of all target mRNAs to provide strong evidence for the further functional verification of these circRNAs. First, GO bioinformatic analysis covers three domains: biological process, cellular component and molecular function. For each part, we showed the following classifications for significantly enriched terms: top 10 counts, top 10-fold enrichment value, top 10 enrichment score value and gene ratio values. The results indicated the following: (1) in terms of biological processes, the most significant enrichment and meaningful terms are metabolic process and RNA splicing,     describing a series of biological events (Figures 5A-D).
(2) In terms of cellular component, the most significant enrichment and meaningful terms are intracellular parts and complexes, describing the components of a cell (Figures 6A-D).
(3) In terms of molecular function, the most significant enrichment and meaningful terms are binding events, suggesting a functional role at the molecular level (Figures 7A-D). Second, the top 10 enrichment score value of significantly enriched pathways based on KEGG pathway analysis is presented in Figure 8A. In addition, the dot plot showed the gene ratio value of the top 10 most significantly enriched pathways ( Figure 8B). The results suggest that these target genes may play a vital role in multiple metabolic activities, and the top two significant pathways are the AMPK signaling pathway and the peroxisome related pathway (Figures 8C,D).

DISCUSSION
In the present study, we found that an array of circRNAs were differentially expressed in the early stage after SCI based on circRNA microarray detection. After validation by qRT-PCR, five candidate circRNAs were selected for further prediction. Next, using bioinformatics tools, circRNA/miRNA/mRNA interactions were constructed according to ceRNA mechanisms, indicating a close relationship between circRNAs, miRNAs, and their target mRNAs. Notably, we predicted that the target mRNAs of these confirmed circRNAs may regulate various biological processes, cellular components and molecular functions, suggesting the functional roles of these circRNAs in pathophysiologic processes. Lastly, target mRNAs were revealed to participate in multiple cell signaling pathways of metabolism, including the AMPK signaling pathway and peroxisome activity, which play a pivotal role in the regulation of cell energy homeostasis, significantly affecting cell survival, death, differentiation, proliferation, and inflammation (Chen et al., 2018b;Di Cara et al., 2018). Taken together, the observations from this study suggest that these altered circRNAs are extensively involved in pathophysiologic processes after traumatic SCI. In reviewing the literature, no studies have focused on screening circRNA expression patterns in any model with traumatic SCI. Given this finding, the present study, to the best of our knowledge, provides the first comprehensive assessment of circRNA expression patterns and subsequent functional prediction after traumatic SCI, which offers strong evidence for further investigation.
Most of the human genome is composed of ncRNAs, which are widely involved in physiological and pathological activities and are closely related to many diseases (Dimartino et al., 2018;Sun et al., 2018). miRNAs and lincRNAs, two important types of ncRNAs, have caused interest due to their multiple functions in diverse pathophysiological events (Nicolas, 2017). Prior studies have noted that miRNAs and lincRNAs, as key factors, regulate the molecular mechanisms behind secondary injury after SCI, including inflammation, blood-brain damage, apoptosis, autophagy, oxidative stress, edema, endoplasmic reticulum stress, and demyelination (Shi et al., 2017). Interestingly, normalizing the level of dysregulated miRNAs and lincRNAs can alleviate the pathological changes after SCI (Hu et al., 2013). circRNAs, a recent rising star of ncRNAs, are endogenously expressed as single-stranded and covalently closed circular molecules that are enriched in mammals. Compared to miRNAs and lincRNAs, they are more stable due to their unique covalently closed loop and specific tertiary structures, which offer more possibilities to act as ideal biomarkers or novel therapeutic targets. Recent growing evidence has shown that circRNAs are implicated in various biological processes and human illness, such as cancers, cardiovascular diseases, and neurological disorders (Floris et al., 2017;Chen et al., 2018a;Hu et al., 2018). As an example, in Alzheimer's disease, Cdr1as is downregulated, thus regulating the downstream target gene and ultimately affecting the progression of disease (Zhao et al., 2016). In terms of trauma to the central nervous system, an report stated for the first time that a total of 192 circRNAs were observed to be differentially expressed (FC ≥ 1.5 and P < 0.05) after traumatic brain injury, 98 of which were upregulated and 94 of which were downregulated, indicating the potential roles of these altered circRNAs in pathophysiologic processes after traumatic brain injury 19 . Another study  provided evidence concerning circRNA expression alteration in exosomes from the brain extracellular space after traumatic brain injury in mice, which broadens the horizon of research on circRNAs. However, no previous studies have profiled circRNA expression in the spinal cord after SCI. As a result, we put forward a hypothesis that circRNAs expression is significantly altered after traumatic SCI in rats. Finally, after validation, the striking results in our study indicated that circRNAs were altered significantly in rat spinal cord after SCI compared to sham controls.
The relationship between circRNAs and their potential target genes was examined to investigate the functions of these dysregulated circRNAs after screening their expression patterns. Mounting evidence has reported that circRNA can serve as ceRNAs for miRNAs, namely, as miRNA sponges (Liu et al., 2017b;Bai et al., 2018). A circRNA may contain multiple miRNA binding sites and may have adsorptive and suppressive effects on miRNAs. Thus, circRNAs can remove the inhibitory effect of miRNAs on their target mRNAs, indirectly upregulating these target mRNAs. As an example, in the plasma of acute ischemic stroke patients and in a mouse stroke model, circDLGAP4 was significantly downregulated, functioning as an endogenous miR-143 sponge to inhibit miR-143 activity and leading to the inhibition of target gene expression (Bai et al., 2018). In addition, in bladder carcinoma, it has been demonstrated that over-expression of circTCF25 could sponge miR-103a-3p and miR-107, increase CDK6 expression, and promote proliferation and migration in vitro and in vivo . Interestingly, several studies stated that a protein can be translated from circRNAs in human cells driven by N6-methyladenosine22, suggesting a possible translatable function of circRNAs (Du et al., 2017;Yang et al., 2017). However, the number of relevant studies is quite limited, and further studies should be conducted to validate these novel results. As a consequence, based on the function of circRNAs as ceRNAs in regulating the activity of corresponding linear mRNAs by binding miRNAs, we first predicted the potential sponging miRNAs for each differentially expressed circRNA by conserved seed sequence matches. We ranked miRNA candidates based on the mirSVR and listed the five highest-ranking miRNA candidates for each circRNA, including the five candidate circRNAs. Then, target-binding mRNAs for miRNA candidates were calculated and filtered with bioinformatics tools. In our study, a total of 5 circRNAs, 60 miRNAs, and 253 mRNAs were included to construct circRNA/mRNA of miRNA interactions, which presented a large interaction network for bioinformatic analysis. These results together lay the groundwork for future research into the specific circRNA/miRNA/mRNA network, which is beneficial for investigating the role of circRNAs in regulating the expression of target genes. We will next validate the selected interaction.
In addition, GO analysis and KEGG pathway analysis were performed to functionally annotate the predicted target mRNAs of circRNAs-miRNAs axes. First, the goal of the Gene Ontology Consortium is to produce a dynamic, controlled vocabulary that can be applied to all eukaryotes, even as knowledge of gene and protein roles in cells is accumulating and changing (Ashburner et al., 2000); thus, three independent ontologies accessible on the World Wide Web 4 are being constructed: biological process, molecular function, and cellular component. Briefly, we found that the target mRNAs are involved in multiple biological processes, cellular signaling pathways, protein activities and gene splicing in rat with SCI. As an example, the carbohydrate metabolic process was found to be one of the most significantly enriched and meaningful terms of biological processes after GO analysis. In addition, one previous study confirmed the central adiposity associations to carbohydrate and lipid metabolism in individuals with complete motor SCI (Gorgey et al., 2011). Thus, the results together suggested the functional roles of the circRNA/miRNA/mRNA network in regulating pathophysiology after SCI. Second, we used KEGG pathway analysis to reveal the potential roles of the target mRNAs regulated by circRNAs-miRNAs axes in diverse biological pathways. According to our annotation, the top two signaling pathways affected by the circRNAs-miRNAs axes were the AMPK signaling pathway and the peroxisome related pathway. Previous studies have explored the relationships between these pathways and pathophysiological processes after SCI. As an example, AMPK is a pivotal regulator of energy homeostasis, and after assessing the influence of longstanding and recent SCI on protein abundance of AMPK isoforms in human skeletal muscle, researchers found that physical/neuromuscular activity is an important determinant of isoform abundance of AMPK (Kostovski et al., 2013). Additionally, another report described the pharmacological activation of peroxisome proliferator-activated receptors due to their anti-inflammatory/antioxidant/anti-excitotoxic/proenergetic profile in SCI (Esposito and Cuzzocrea, 2011). After pathway analysis, we performed pathway diagrams of the AMPK signaling pathway and the peroxisome related pathway (Figures 8A,B), presenting the most significant core downstream pathways of the circRNAs-miRNAs axes directly.
However, several questions remain to be answered. First and most importantly, the reader should bear in mind that this study was unable to validate the circRNA/miRNA/mRNA network predicted by bioinformatics, and our future work is to identify their functions. Next, the relatively small group size for the PCR experiments may limit the interpretation for candidate circRNAs and subsequent bioinformatics. Thus, the enlarged group size is needed to improve the scientific rigor. Then, it was beyond the scope of this study to examine the potential effect of different time points on the possible dynamic changes in circRNA expression patterns in the contused spinal cords of rats with SCI. Last but not least, we strongly recommend that future research should focus on different animal models and types of SCI to study the expression and functions of circRNAs.
Taken together, the current findings of this study were greatly significant in at least two major respects. On one hand, this study explored, for the first time, the significant alteration of circRNA expression profiles after traumatic SCI in rats. On the other hand, these differentially expressed circRNAs were predicted to closely correlate with post-SCI pathophysiology processes based on bioinformatics analysis. Despite its exploratory nature, this study offered deep insights into many possible treatment targets of SCI by regulating circRNAs.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.