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

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 (Fan et al., 2017), 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 (Liu et al., 2017a). 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.

Materials and Methods

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 (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¹). 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²) 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³). 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.

Results

Expression Pattern of CircRNAs in Spinal Cord Tissues After SCI

Prior to injury, the BBB scores of all rats were 21. Rats in the sham group remained at 21, indicating the integrity of the spinal cord during the experimental period. However, in the SCI group, no hindlimb motor performance was observed immediately after 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.

TABLE 1

Table 1. The top 10 up-regulated and down-regulated circRNAs ranked by fold changes after SCI.

FIGURE 1

Figure 1. Hierarchical clustering of circRNA expression. C101/103/104 refers to the cord samples from the sham control group, and T113/114/116 refers to the cord samples from the SCI group. (A) Hierarchical clustering analysis included all 13279 circRNAs between the sham control and the SCI groups. (B) Hierarchical clustering analysis included differentially expressed circRNAs (fold change ≥ 2; P < 0.05) between the sham control and the SCI groups.

FIGURE 2

Figure 2. Differences in the circRNA expression profiles between the two groups. (A) The scatter plot showed the differences in circRNA expression between the SCI and sham groups. The values of the X and Y axes in the scatter plot are the normalized signal values of the samples (log2 scaled) or the averaged normalized signal values of groups of samples (log2 scaled). The green lines are fold change lines. The circRNAs above the top green line and below the bottom green line indicated more than twofold changes of circRNAs between the two compared samples. (B) Volcano plots show the differentially expressed circRNAs with statistical significance (fold change ≥ 2; P < 0.05). The vertical lines correspond to 2.0-fold up and down, respectively, and the horizontal line represents a P of 0.05; the red point in the plot represents differentially expressed circRNAs with statistical significance. (C) The distribution of differentially expressed circRNAs in chromosomes is presented, showing that the dysregulated circRNAs stem from every chromosome. (D) The bar diagram of circRNA categories based on gene sources is shown, revealing that most of the circRNAs altered after SCI are exonic.

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).

TABLE 2

Table 2. The five highest-ranking miRNA candidates for top 10 up-regulated and down-regulated 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

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.

TABLE 3

Table 3. Comparison for candidate circRNAs expression in microarray and PCR.

TABLE 4

Table 4. Sequences of primers used for qRT-PCR assay.

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 can serve as ceRNA for miRNAs, and miRNAs usually inhibit target mRNAs; therefore, circRNAs may indirectly upregulate target mRNAs by inhibiting the negative regulation of miRNAs.

FIGURE 4

Figure 4. The circRNA/miRNA/mRNA network analysis. The network included the 5 circRNAs, 60 miRNAs, and 253 mRNAs (Nodes with red color are miRNAs; nodes with light-blue color are mRNAs; nodes with brown color are 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).

FIGURE 5

Figure 5. The GO annotations for biological process of target mRNAs regulated by the five candidate circRNAs. (A) The pie chart above shows the top 10 counts of the significant enrichment terms. (B) The bar plot shows the top 10 enrichment score values of the significant enrichment terms. (C) The bar plot shows the top 10-fold enrichment values of the significant enrichment terms. (D) The dot plot shows the gene ratio values of the top 10 most significant enrichment terms. Gene ratio: the GOID’s gene ratio value, which equals (Count/List.Total).

FIGURE 6

Figure 6. The GO annotations for cellular component of target mRNAs regulated by the five candidate circRNAs. (A) The pie chart above shows the top 10 counts of the significant enrichment terms. (B) The bar plot shows the top 10 Enrichment Score values of the significant enrichment terms. (C) The bar plot shows the top 10-fold enrichment values of the significant enrichment terms. (D) The dot plot shows the gene ratio values of the top 10 most significant enrichment terms. Gene ratio: GOID’s gene ratio value, which equals (Count/List.Total).

FIGURE 7

Figure 7. The GO annotations for molecular function of target mRNAs regulated by the five candidate circRNAs. (A) The pie chart above shows the top 10 counts of the significant enrichment terms. (B) The bar plot shows the top 10 enrichment score values of the significant enrichment terms. (C) The bar plot shows the top 10-fold enrichment values of the significant enrichment terms. (D) The dot plot shows the gene ratio values of the top 10 most significant enrichment terms. Gene ratio: the GOID’s gene ratio value, which equals (Count/List.Total).

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).

FIGURE 8

Figure 8. KEGG pathway analysis of target mRNAs regulated by the five candidate circRNAs. (A) The bar plot shows the top 10 enrichment score values of the significantly enriched pathway. (B) The dot plot shows the gene ratio values of the top 10 most significantly enriched pathways. (C) The AMP-activated protein kinase (AMPK) signaling pathway and (D) the peroxisome related pathway are shown. Yellow marked nodes are associated with downregulated genes, orange marked nodes are associated with upregulated or only whole dataset genes, and green nodes have no significance.

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¹⁹. Another study (Zhao et al., 2018) 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 (Zhong et al., 2016). 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⁴ 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/pro-energetic 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.

Author Contributions

CQ, C-BL, and M-LY performed the research. J-JL designed the research study. D-GY and FG contributed essential reagents or tools. CZ and L-JD analyzed the data. CQ wrote the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by a grant from the Special Fund for Basic Scientific Research of Central Public Research Institutes (2018CZ-1), the National Natural Science Foundation of China (General Program 81870979 and 81271366), Basic Scientific Research Foundation of China Rehabilitation Research Center (2018ZX-30), the Scientific Research Foundation of CRRC (2012C-1), the Beijing Science and Technology Major Project (D161100002816004), and the Capital Health Research and Development of Special (2018-1-6011).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

Special thanks to Pei-Pei Xu for her animal care.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnmol.2018.00497/full#supplementary-material

Abbreviations

AMPK, AMP-activated protein kinase; ceRNAs, competing endogenous RNAs; circRNAs, circular RNAs; FCs, fold changes; FDR, false discovery rate; GAPDH, glyceraldehyde 3-phosphate dehydrogenase; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes database; lncRNA, long non-coding RNA; miRNA, microRNA; MREs, miRNA response elements; ncRNAs, non-coding RNAs; qRT-PCR, quantitative reverse transcription-PCR; SCI, spinal cord injury.

Footnotes

1. ^https://sourceforge.net/projects/mev-tm4/files/mev-tm4/
2. ^http://www.geneontology.org
3. ^www.genome.jp/kegg
4. ^http://www.geneontology.org

References

Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene ontology: tool for the unification of biology. the gene ontology consortium. Nat. Genet. 25, 25–29. doi: 10.1038/75556

PubMed Abstract | CrossRef Full Text | Google Scholar

Bai, Y., Zhang, Y., Han, B., Yang, L., Chen, X., Huang, R., et al. (2018). Circular RNA DLGAP4 ameliorates ischemic stroke outcomes by targeting mir-143 to regulate endothelial-mesenchymal transition associated with blood-brain barrier integrity. J. Neurosci. 38, 32–50. doi: 10.1523/JNEUROSCI.1348-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Basso, D. M., Beattie, M. S., and Bresnahan, J. C. (1996). Graded histological and locomotor outcomes after spinal cord contusion using the NYU weight-drop device versus transection. Exp. Neurol. 139, 244–256. doi: 10.1006/exnr.1996.0098

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Zhang, L., Su, Y., and Zhang, X. (2018a). Screening potential biomarkers for colorectal cancer based on circular RNA chips. Oncol. Rep. 39, 2499–2512. doi: 10.3892/or.2018.6372

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X. L., Wang, Y., Peng, W. W., Zheng, Y. J., Zhang, T. N., Wang, P. J., et al. (2018b). Effects of interleukin-6 and IL-6/AMPK signaling pathway on mitochondrial biogenesis and astrocytes viability under experimental septic condition. Int. Immunopharmacol. 59, 287–294. doi: 10.1016/j.intimp.2018.04.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Di Cara, F., Sheshachalam, A., Braverman, N. E., Rachubinski, R. A., and Simmonds, A. J. (2018). Peroxisome-mediated metabolism is required for immune response to microbial infection. Immunity 48, 832–833. doi: 10.1016/j.immuni.2018.03.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Dimartino, D., Colantoni, A., Ballarino, M., Martone, J., Mariani, D., Danner, J., et al. (2018). The Long non-coding RNA lnc-31 interacts with rock1 MRNA and mediates its YB-1-dependent translation. Cell Rep. 23, 733–740. doi: 10.1016/j.celrep.2018.03.101

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, W. W., Zhang, C., Yang, W., Yong, T., Awan, F. M., and Yang, B. B. (2017). Identifying and characterizing circRNA-protein interaction. Theranostics 7, 4183–4191. doi: 10.7150/thno.21299

PubMed Abstract | CrossRef Full Text | Google Scholar

Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA targets in drosophila. Genome Biol. 5:R1.

Google Scholar

Esposito, E., and Cuzzocrea, S. (2011). Targeting the peroxisome proliferator-activated receptors (PPARs) in spinal cord injury. Expert Opin. Ther. Targets 15, 943–959. doi: 10.1517/14728222.2011.581231

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, X., Wang, J. Z., Lin, X. M., and Zhang, L. (2017). Stem cell transplantation for spinal cord injury: a meta-analysis of treatment effectiveness and safety. Neural Regen. Res. 12, 815–825. doi: 10.4103/1673-5374.206653

PubMed Abstract | CrossRef Full Text | Google Scholar

Floris, G., Zhang, L., Follesa, P., and Sun, T. (2017). Regulatory role of circular RNAs and neurological disorders. Mol. Neurobiol. 54, 5156–5165. doi: 10.1007/s12035-016-0055-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman, R. C., Farh, K. K., Burge, C. B., and Bartel, D. P. (2009). Most mammalian mRNAs are conserved targets of microRNAs. Genome Res. 19, 92–105. doi: 10.1101/gr.082701.108

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorgey, A. S., Mather, K. J., and Gater, D. R. (2011). Central adiposity associations to carbohydrate and lipid metabolism in individuals with complete motor spinal cord injury. Metabolism 60, 843–851. doi: 10.1016/j.metabol.2010.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Hextrum, S., and Bennett, S. (2018). A critical examination of subgroup analyses: the national acute spinal cord injury studies and beyond. Front. Neurol. 9:11. doi: 10.3389/fneur.2018.00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, A. M., Li, J. J., Sun, W., Yang, D. G., Yang, M. L., Du, L. J., et al. (2015). Myelotomy reduces spinal cord edema and inhibits aquaporin-4 and aquaporin-9 expression in rats with spinal cord injury. Spinal Cord 53, 98–102. doi: 10.1038/sc.2014.209

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, J. Z., Huang, J. H., Zeng, L., Wang, G., Cao, M., and Lu, H. B. (2013). Anti-apoptotic effect of microRNA-21 after contusion spinal cord injury in rats. J. Neurotrauma 30, 1349–1360. doi: 10.1089/neu.2012.2748

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, W., Bi, Z. Y., Chen, Z. L., Liu, C., Li, L. L., Zhang, F., et al. (2018). Emerging landscape of circular RNAs in lung cancer. Cancer Lett. 427, 18–27. doi: 10.1016/j.canlet.2018.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, M., Zhong, Z., Lv, M., Shu, J., Tian, Q., and Chen, J. (2016). Comprehensive analysis of differentially expressed profiles of lncRNAs and circRNAs with associated co-expression and ceRNA networks in bladder carcinoma. Oncotarget 7, 47186–47200. doi: 10.18632/oncotarget.9706

PubMed Abstract | CrossRef Full Text | Google Scholar

Hwang, T., Park, C. K., Leung, A. K., Gao, Y., Hyde, T. M., Kleinman, J. E., et al. (2016). Dynamic regulation of RNA editing in human brain development and disease. Nat. Neurosci. 19, 1093–1099. doi: 10.1038/nn.4337

PubMed Abstract | CrossRef Full Text | Google Scholar

Khuyagbaatar, B., Kim, K., and Kim, Y. H. (2015). Conversion equation between the drop height in the New York university impactor and the impact force in the infinite horizon impactor in the contusion spinal cord injury model. J. Neurotrauma 32, 1987–1993. doi: 10.1089/neu.2015.3875

PubMed Abstract | CrossRef Full Text | Google Scholar

Kostovski, E., Boon, H., Hjeltnes, N., Lundell, L. S., Ahlsen, M., Chibalin, A. V., et al. (2013). Altered content of AMP-activated protein kinase isoforms in skeletal muscle from spinal cord injured subjects. Am. J. Physiol. Endocrinol. Metab. 305, E1071–E1080. doi: 10.1152/ajpendo.00132.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Larionov, A., Krause, A., and Miller, W. (2005). A standard curve based method for relative real time PCR data processing. BMC Bioinformatics 6:62. doi: 10.1186/1471-2105-6-62

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, C., Zhang, C., Yang, J., Geng, X., Du, H., Ji, X., et al. (2017a). Screening circular RNA expression patterns following focal cerebral ischemia in mice. Oncotarget 8, 86535–86547. doi: 10.18632/oncotarget.21238

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Liu, T., Wang, X., and He, A. (2017b). Circles reshaping the RNA world: from waste to treasure. Mol. Cancer 16:58. doi: 10.1186/s12943-017-0630-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, H. R. (2017). lncRNA-Map2k4 sequesters miR-199a to promote FGF1 expression and spinal cord neuron growth. Biochem. Biophys. Res. Commun. 490, 948–954. doi: 10.1016/j.bbrc.2017.06.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Manohar, A., Foffani, G., Ganzer, P. D., Bethea, J. R., and Moxon, K. A. (2017). Cortex-dependent recovery of unassisted hindlimb locomotion after complete spinal cord injury in adult rats. eLife 6:e23532. doi: 10.7554/eLife.23532

PubMed Abstract | CrossRef Full Text | Google Scholar

Memczak, S., Jens, M., Elefsinioti, A., Torti, F., Krueger, J., Rybak, A., et al. (2013). Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 495, 333–338. doi: 10.1038/nature11928

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagoshi, N., and Okano, H. (2017). Applications of induced pluripotent stem cell technologies in spinal cord injury. J. Neurochem. 141, 848–860. doi: 10.1111/jnc.13986

PubMed Abstract | CrossRef Full Text | Google Scholar

Nicolas, F. E. (2017). Role of ncRNAs in development, diagnosis and treatment of human cancer. Recent Pat. Anticancer Drug Discov. 12, 128–135. doi: 10.2174/1574892812666170105113415

PubMed Abstract | CrossRef Full Text | Google Scholar

Pang, X., Tang, Y., and Zhang, D. (2016). Role of miR-145 in chronic constriction injury in rats. Exp. Ther. Med. 12, 4121–4127. doi: 10.3892/etm.2016.3900

PubMed Abstract | CrossRef Full Text | Google Scholar

Pasquinelli, A. E. (2012). MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship. Nat. Rev. Genet. 13, 271–282. doi: 10.1038/nrg3162

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, S., Yang, X., Li, X., Wang, J., Gao, Y., Shang, R., et al. (2015). Circular RNA: a new star of noncoding RNAs. Cancer Lett. 365, 141–148. doi: 10.1016/j.canlet.2015.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Salgado-Somoza, A., Zhang, L., Vausort, M., and Devaux, Y. (2017). The circular RNA MICRA for risk stratification after myocardial infarction. Int. J. Cardiol. Heart Vasc. 17, 33–36. doi: 10.1016/j.ijcha.2017.11.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell 146, 353–358. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheff, S. W., Saucier, D. A., and Cain, M. E. (2002). A statistical method for analyzing rating scale data: the BBB locomotor score. J. Neurotrauma 19, 1251–1260. doi: 10.1089/08977150260338038

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, Z., Zhou, H., Lu, L., Li, X., Fu, Z., Liu, J., et al. (2017). The roles of microRNAs in spinal cord injury. Int. J. Neurosci. 127, 1104–1115. doi: 10.1080/00207454.2017.1323208

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, N., Zhang, G., and Liu, Y. (2018). Long non-coding RNA XIST sponges miR-34a to promotes colon cancer progression via Wnt/beta-catenin signaling pathway. Gene 665, 141–148. doi: 10.1016/j.gene.2018.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, B., Wang, Y., Lin, Y., Zhao, C. C., Mao, Q., Feng, J. F., et al. (2018). Circular RNA Expression profiles alter significantly after traumatic brain injury in rats. J. Neurotrauma 35, 1659–1666. doi: 10.1089/neu.2017.5468

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, D. G., Li, J. J., Gu, R., Yang, M. L., Zhang, X., Du, L. J., et al. (2013). Optimal time window of myelotomy in rats with acute traumatic spinal cord injury: a preliminary study. Spinal Cord 51, 673–678. doi: 10.1038/sc.2013.56

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, M. L., Li, J. J., Gao, F., Du, L. J., Zhao, H. P., Wang, Y. M., et al. (2014). A preliminary evaluation of the surgery to reconstruct thoracic breathing in patients with high cervical spinal cord injury. Spinal Cord 52, 564–569. doi: 10.1038/sc.2014.64

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Fan, X., Mao, M., Song, X., Wu, P., Zhang, Y., et al. (2017). Extensive translation of circular RNAs driven by N(6)-methyladenosine. Cell Res. 27, 626–641. doi: 10.1038/cr.2017.31

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, R., Zhou, J., Dong, X., Bi, C., Jiang, R., Dong, J., et al. (2018). Circular RNA expression alteration in exosomes from the brain extracellular space after traumatic brain injury in mice. J. Neurotrauma 35, 2056–2066. doi: 10.1089/neu.2017.5502

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Alexandrov, P. N., Jaber, V., and Lukiw, W. J. (2016). Deficiency in the ubiquitin conjugating enzyme UBE2A in alzheimer’s disease (AD) is Linked to deficits in a natural circular miRNA-7 sponge (circRNA; ciRS-7). Genes 7:116. doi: 10.3390/genes7120116

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, Z., Lv, M., and Chen, J. (2016). Screening differential circular RNA expression profiles reveals the regulatory role of circTCF25-miR-103a-3p/miR-107-CDK6 pathway in bladder carcinoma. Sci. Rep. 6:30919. doi: 10.1038/srep30919

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, H. J., Wang, L. Q., Wang, D. B., Yu, J. B., Zhu, Y., Xu, Q. S., et al. (2018a). Long non-coding RNA MALAT1 contributes to inflammatory response of microglia following spinal cord injury via modulating miR-199b/IKKbeta/NF-kappaB signaling pathway. Am. J. Physiol. Cell Physiol. 315, C52–C61. doi: 10.1152/ajpcell.00278.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, R., Wu, Y., Wang, W., Su, W., Liu, Y., Wang, Y., et al. (2018b). Circular RNAs (circRNAs) in cancer. Cancer Lett. 425, 134–142. doi: 10.1016/j.canlet.2018.03.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, H., Xie, R., Liu, X., Shou, J., Gu, W., Gu, S., et al. (2017a). MicroRNA-494 improves functional recovery and inhibits apoptosis by modulating PTEN/AKT/mTOR pathway in rats after spinal cord injury. Biomed. Pharmacother. 92, 879–887. doi: 10.1016/j.biopha.2017.05.143

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, L. L., Li, M. Q., He, F., Zhou, S. B., and Jiang, W. (2017b). Mitochondria targeted peptide attenuates mitochondrial dysfunction, controls inflammation and protects against spinal cord injury-induced lung injury. Cell Physiol. Biochem. 44, 388–400. doi: 10.1159/000484919

PubMed Abstract | CrossRef Full Text | Google Scholar