Young adult C57BL/6 mice (female, 25–30 g, 8 weeks) were housed for at least seven days in a controlled temperature environment (23°C) on a 12 h light/dark cycle before the operation. The mice were allowed ad libitum access to water and food except on the day before the operation, on which the mice were allowed only ad libitum access to water. All animal procedures were approved by the Institutional Animal Care and Use Committee of the Sir Run Run Shaw Hospital of Zhejiang University School of Medicine, People's Republic of China.

The SCI model used Allen's method, as described in our previous study (15). Twelve mice were anesthetized intraperitoneally by 1% pentobarbital sodium. They were randomly divided into two groups: SCI group (n = 6) and sham group (n = 6). The dorsal surface of thoracic spinal cord segments was exposed, taking care to avoid any dural tears at the T9-T10 level after laminectomy. Subsequently, a New York University (NYC) II struck instrument with a 1-mm-diameter impactor from a height of 25 mm was used to induce a moderate contusion injury of the spinal cord at T10. The wound was then closed by suturing the muscle layer using nylon sutures. The same procedure without the contusion injury was performed in the sham group.

Three days after SCI, six mice (SCI group, n = 3; sham group, n = 3) were anesthet ized. One cm long segments at the T9-T10 spinal cord enclosing the injured site were rapidly frozen in liquid nitrogen and transported to Shanghai Biotechnology Corporation for microarray analysis on dry ice. Total RNA was extracted and purified using the mirVana™ miRNA Isolation Kit (Cat#AM1561, Ambion, Austin, TX, US) following the manufacturer's instructions and the RIN was assessed on the Agilent Bioanalyzer 2,100 (Agilent technologies, Santa Clara, CA, US).

Total RNA was amplified and labeled using the Low Input Quick Amp Labeling Kit, One-Color (Cat.#5190-2305, Agilent technologies, Santa Clara, CA, US) following the manufacturer's instructions. Labeled cRNAs were purified using the RNeasy mini kit (Cat. #74106, QIAGEN, GmBH, Germany). Each slide was hybridized with 1.65 μg of Cy3-labeled cRNA using the Gene Expression Hybridization Kit (Cat. #5188-5242, Agilent technologies, Santa Clara, CA, US) in a Hybridization Oven (Cat. #G2545A, Agilent technologies, Santa Clara, CA, US). Slides were then washed in staining dishes according to the manufacturer's instructions (Cat. #121, Thermo Shandon, Waltham, MA, US) with the Gene Expression Wash Buffer Kit (Cat. #5188-5327, Agilent technologies, Santa Clara, CA, US) after 17 h hybridization. Slides were scanned on a Agilent Microarray Scanner (Cat#G2565CA, Agilent technologies, Santa Clara, CA, US) with default settings, as follows: Dye channel: Green, scan resolution = 3μm, PMT 100%, 20 bit. Raw data were normalized using the quantile algorithm of the limma package in R. Data were extracted using the Feature Extraction software 12.0 (Agilent technologies, Santa Clara, CA, US).

Scanned images were imported into the Agilent Feature Extraction software for raw data extraction. Using the R software package limma, quantile normalization of raw data and data processing was conducted and low intensity probe filtering was performed. The log2-ratio was used for quantile normalization. Via fold change filtering or volcano plot filtering, differentially expressed circRNAs with statistical significance were identified. Hierarchical clustering was performed to highlight differential circRNA expression patterns among different samples. Significantly differentially expressed circRNAs were selected as follows: Fold change ≥2, p-value <0.05. We entered the host genes of differentially expressed circRNAs into the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov) for annotation and functional analysis. GO analysis including the categories of biological process, cellular component, and molecular function was performed on account of GO (http://www.geneontology.org). The significance of GO term enrichment was indicated by fold change and P-value (FC ≥2 and P < 0.05 were considered statistically significant). We also used pathway analysis to classify differentially expressed genes on account of the Kyoto Encyclopedia of Genes and Genomes (KEGG) database to predict pathways of interest. mRNAs related to vascular endothelial proliferation, migration, and angiogenesis from the mRNAs which changed more than 2-fold and a P-value < 0.05 were chosen. We then assessed circRNA correlation with angiogenesis via cis-acting mRNA networks.

GO analysis, which was performed to explore the potential functions of parental genes in terms of the circRNAs, covered three different aspects: Biological process (BP), cellular component (CC), and molecular function (MF). The top 10 enriched GO terms and pathways among the two groups were ranked by enrichment score [-log10 (P-value)] as identified by the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://www.david.abcc.ncifcrf.gov/).

qRT-PCR was used to validate circRNA expression patterns following microarray detection. Briefly, using the ReverTra Ace qPCR Kit (TOYOBO, FSQ-101) and High Capacity cDNA Reverse Transcription Kit with RNase Inhibitor (ABI, 4374966), we reverse-transcribed total RNA. qRT-PCR was conducted in the 7,900 HT Sequence Detection System (ABI, USA) and Quant Studio 5 Real-Time PCR System (ABI, USA) using Power SYBR Green PCR Master Mix (ABI, 4368708). We used Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as an internal control to normalize the data. The primers of these circRNAs and GAPDH are presented in Table 1.

circRNA microarray expression profiling and data bioinformatics analysis were performed by Shanghai Biotechnology Corporation, Shanghai, People's Republic of China. Based on the expression signals, we built a circRNA-miRNA-mRNA interaction network using a regression model and seed sequence matching analysis. The circRNA-mRNA interactions were predicted using the principle of high related coefficient and at least one perfect seed-matching sequence.

Three days after SCI, six mice (SCI group, n = 3; sham group, n = 3) were euthanized via an overdose of 1% pentobarbital sodium and the spinal cord lesion site was extracted and flash-frozen at −80°C. The spinal cords were triturated using Tissue Prep (Gering Scientific Instruments, Beijing, China) and cells were dissolved in WB and IP lysis buffer containing phenylmethanesulfonylfluoride (PMSF) and phosphatase inhibitor and centrifuged at 1.3 × 10⁴ g for 15 min at 4°C. Subsequently, the protein concentration was quantified using a BCA kit (Beyotime Institute of Biotechnology, Jiangsu, China). All samples (20 μl) were separated on a SurePAGE^TM (Cat.M00653, GenScript®, USA) and transferred onto nitrocellulose membranes. Membranes were blocked in 5% skimmed milk powder dissolved in TBS-T for 1 h. Mouse anti-CD31 (Cat# 12242S, ABclonal, 1:1000), rabbit anti-VEGF (Cat# 3495S, Cell Signaling Technology, 1:1000), and rabbit anti-VEGFR2 (Cat# 4108S, ABclonal, 1:1000) were used for probing interested proteins at 4°C overnight. Membranes were incubated with infrared-labeled secondary antibodies (Li-COR Biosciences) after several washes with TBS-T. In order to visualize the immunoblot bands, an Odyssey infrared imaging system (LI-COR® Biosciences, NE, USA) was used. Band intensity was analyzed with an Image Studio Ver 5.2 system and standardized to the β-Actin internal standard.

All quantitative data are presented as mean ± standard deviation (mean ± SD). Student's t-test was used to calculate the statistical significance between the groups. All statistical analyses were performed using the SPSS 20.0 software (SPSS, Chicago, IL). The correlation of two variables was assessed using the Pearson product-moment correlation coefficient. Genes changing equal to or more than 2-fold with P < 0.05 were considered statistically significant.

The microarrays used in the study contained probes for 37,852 mouse circRNAs, of which 9364 circRNAs were found on circBase, 378 circRNAs were found on Deepbase, and the remaining were found on Shbio. As shown in the box plot in Figure 1A, the distribution of circRNA expression profiles was not different in samples, indicating that post-injury changes observed in the levels of individual circRNAs were not due to chance. Differentially expressed circRNAs were filtered by fold change and p-value (>2-fold change, p < 0.05), which are illustrated as red and blue dots on the scatter plot in Figure 1B and volcano plot in Figure 1C. Hierarchical clustering was performed to highlight distinguishable circRNA expression patterns between the two groups (Figure 1D). The altered circRNAs were clearly separated into two clusters, indicating that the samples had good intra-group consistency, while circRNA expression levels in the SCI group were significantly different compared to the sham group. Compared with the sham group, 1288 circRNAs were significantly altered (>2-fold change, p < 0.05) in the spinal cord after SCI, of which 991 were upregulated and 297 were downregulated. The top 20 altered circRNAs are listed in Table 1. The distributions of altered circRNAs in mice revealed that the upregulated circRNAs were transcribed from all chromosomes, while downregulated circRNAs were transcribed from all chromosomes except chrY (Figure 2A).

In order to verify the microarray data, real-time PCR was used to evaluate the expression of two randomly selected circRNAs. Primers used for qRT-PCR of these circRNAs are listed in Table 2. Compared with the sham group, circRNA.7079 was upregulated in the injured group (Figure 2B), while circRNA.12884 was downregulated significantly after SCI (Figure 2C).

In order to analyze the pathways associated with endothelial proliferation, migration, and angiogenesis and assess the numbers of related circRNAs, we performed KEGG analyses. In order to reveal the mRNAs and circRNAs correlated with vascular endothelial proliferation, migration, and angiogenesis, we performed GO analyses. We conducted GO and KEGG pathway analysis for the differentially expressed circRNAs based on their host genes in order to investigate potential biological functions. The GO analysis indicated that “biological regulation,” “cellular process,” and “single-organism process response” were the main related biological processes involved in SCI (Figure 3A). In terms of cellular components, “cell,” “cell part,” and “organelle” were highly activated. Most of the molecular functions activated after SCI were “binding,” “catalytic activity,” and “molecular function regulator.” KEGG analysis showed that the top 30 pathways affected by the altered circRNAs might be involved in the progression of SCI (Figure 3B). The host-genes of circRNAs were predominantly enriched for glycosaminoglycan biosynthesis and extracellular matrix (ECM)-receptor interaction pathways. The most significant pathway was the ECM-receptor interaction, which also had a high enrichment score, reflecting that this pathway was affected in the spinal cord after SCI more than in other diseases or conditions.

The competing-endogenous RNA (ceRNA) network analysis proposes that RNA transcripts, both coding and non-coding, crosstalk with and coregulate each other using microRNA response elements (MREs). The ceRNA network analysis tremendously expands functional information of coding and non-coding RNAs. Mounting evidence has shown that various types of RNAs, including pseudogenes, long non-coding RNAs, circular RNAs, and messenger RNAs, can function as ceRNAs in distinct physiological and pathophysiological states. We screened the top 20 differently expressed circRNAs after SCI and identified three of them for which related mRNAs could be predicted via ceRNA analysis, circRNA.7079, circRNA.7078, and circRNA.6777. The integrated circRNA-miRNA-mRNA network is shown in Figure 4.

We were interested in the expression profiles of circRNA which correlated with vascular endothelial proliferation, migration, and angiogenesis after SCI, thus datasets concerning other functions were excluded. Firstly, to confirm angiogenesis-related protein expression in the spinal cord, western blotting was performed to detect the expression of angiogenesis-related proteins CD31, VEGF, and VEGFR2 3 days after SCI (Figure 5A). The expression levels of CD31, VEGF, and VEGFR2 were strongly decreased at 3 days after SCI (0.390 ± 0.055-fold of sham, 0.255 ± 0.044-fold of sham, and 0.261 ± 0.034-fold of sham, P < 0.001 for all; Figures 5B–D), revealing a reduction of angiogenesis. Next, via KEGG analysis, we analyzed the pathways associated with endothelial proliferation, migration, and angiogenesis and assessed the numbers of related circRNAs (Figure 5E). The host-genes of circRNAs were predominantly enriched for ECM-receptor interaction (n = 97) and the PI3K-Akt signaling pathway (n = 85), as shown in Figure 5E. In the VEGF signaling pathway, several host-genes of circRNAs, such as PLCy, SPK, Ras, COX2, MAPKAFK, and HSP27 were upregulated, while PKC and NEAT were downregulated (KEGG mmu04370, Figure 5F). VEGFR2 is the main receptor of VEGFA which then triggers vascular endothelial proliferation, migration, and angiogenesis through the VEGF signaling pathway, as shown in KEGG mmu04370.

Next, the mRNAs and circRNAs related to vascular endothelial proliferation, migration, and angiogenesis were analyzed using GO analyses. This revealed that 191 mRNAs (exhibiting a >2-fold change and p < 0.01) were correlated with vascular endothelial proliferation, migration, and angiogenesis; such as “vascular endothelial growth factor production (GO:0010573),” “regulation of vascular endothelial growth factor production (GO:0010574),” “positive regulation of angiogenesis (GO:0045766),” “regulation of angiogenesis (GO:0045765),” and “angiogenesis (GO:0001525),” “sprouting angiogenesis (GO:0002040),” amongst other terms (Figures 6A–C). Seventy-five circRNAs (>2-fold change, p < 0.01) were speculated to be correlated with these three biological processes using ceRNA analysis (Figure 6D). In addition, we found 127 mRNAs (>2-fold change, p < 0.01) correlated with vascular endothelial proliferation, migration, and angiogenesis based on KEGG analyses, such as “ECM-receptor interaction (PATHWAYID: mmu04512),” “Ras signaling pathway (PATHWAYID: mmu04014),” “PI3K-Akt signaling pathway (PATHWAYID: mmu04151),” and “VEGF signaling pathway (PATHWAYID: mmu04370),” amongst other terms (Figure 7A). Forty-six circRNAs (>2-fold change, p < 0.01) were predicted to correlate with these mRNAs using ceRNA analysis (Figure 7B). The top 10 upregulated circRNAs were circRNA.7079, circRNA.7078, circRNA.27570, circRNA.26810, circRNA.5566, circRNA.27348, circRNA.6727, circRNA.26807, circRNA.15929, and circRNA.27347. Among these top 10 circRNAs, circRNA.7079, circRNA.7078, circRNA. 27348, and circRNA.15929 were predicted to be involved in vascular endothelial proliferation, migration, and angiogenesis. circRNA.27570 and circRNA.26810 were predicted to be correlated with both vascular endothelial proliferation and migration, while circRNA.5566 and circRNA.6727 were predicted to be involved in vascular endothelial proliferation only. Thus, these circRNAs were predominantly involved in angiogenesis signaling pathways, such as the PI3K-Akt signaling pathway, the STAT3 signaling pathway, the TCR signaling pathway, and the Wnt/β-catenin signaling pathway (Table 3).