Screening Circular RNA Related to Vascular Endothelial Proliferation, Migration and Angiogenesis in Spinal Cord Injury by RNA Sequencing

Background: Traumatic spinal cord injury (SCI) causes high rates of worldwide morbidity because of the complex secondary injury. Circular RNAs (circRNAs) are a novel class of endogenous non-coding RNAs, which have recently been recognized as important regulators of gene expression and pathological processes. In this study, we have attempted to elucidate the expression proles of circRNAs in a mouse model of SCI and comprehensively understand vascular endothelial proliferation, migration and angiogenesis in the early stage of secondary injury. Methods: Deep RNA sequencing (RNA-seq) and bioinformatic analysis including GO enrichment analysis, KEGG pathway analysis and circRNA-miRNA-mRNA network construction were performed to investigate the expression patterns of circRNAs in mouse spinal cord after SCI (n= 3 per group) for three days and explore the differentially expressed circRNAs related to vascular endothelial proliferation, migration and angiogenesis. Results: Total of 1288 circRNAs were altered (>2-fold change, p<0.05) in the spinal cord after SCI, including 991 were upregulated and 297 were downregulated. Meanwhile we constructed a circRNA-mRNA network to predict their functions for circRNAs can act as “miRNA sponges”,. We next analyzed the altered circRNAs related to vascular endothelial proliferation, migration and angiogenesis by GO and KEGG analyses. 121 circRNAs were found to correlating to vascular endothelial proliferation (cid:0) migration and angiogenesis in spinal cord after SCI. Conclusions: Our results reveal that circRNAs locally regulate their related protein-gene expression and play key roles in the vascular endothelial proliferation, migration and angiogenesis of SCI.


Introduction
Traumatic spinal cord injury (SCI) is a devastating condition often affecting young and healthy individuals around the world. This debilitating condition not only creates enormous physical and emotional cost to individuals but also is a signi cant nancial burden to society at large [1]. Patients with SCI often experience devastating neurological 2.2. SCI model Twelve mice were treated with intraperitoneal injection of 1% pentobarbital sodium for deep anesthetization. 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 after laminectomy at the T9-T10 level. Then, a New York University (NYC) 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. After that, by suturing the muscle layer using nylon sutures, the wound was closed. The same procedure without the contusion injury was performed in the sham group.

RNA extraction and puri cation
Three days after SCI, six mice (SCI group, n = 3; sham group, n = 3) were anesthetized and the T9-T10 spinal cord 1 cm long segments centered on and enclosing the injured site were frozen quickly in liquid nitrogen and transported on dry ice to Shanghai Biotechology Corporation for microarray analysis. Following the manufacturer's instructions, total RNA was extracted and puri ed using mirVana™ miRNA Isolation Kit (Cat#AM1561, Ambion, Austin, TX, US) and checked for a RIN number to inspect RNA integration by an Agilent Bioanalyzer 2100 (Agilent technologies, Santa Clara, CA, US).

Bioinformatics analysis
Scanned images were imported into the Agilent Feature Extraction software for raw data extraction. Using the R software package, quantile normalization of raw data and data processing was conducted. subsequently, low intensity ltering was performed. The log2-ratio was used for quantile normalization. Through fold change ltering or volcano plot ltering,differentially expressed circRNAs with statistical signi cance were identi ed. Hierarchical clustering was performed to show the distinguishable circRNA expression pattern among different samples. As signi cantly differentially expressed circRNAs, they were selected(Fold changes ≥ 2, P value < 0.05 was considered statistically signi cant). We put the host genes of differentially expressed circRNAs into the Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov) for the annotation and functional analysis. GO analysis included the categories of biological process, cellular component and molecular function were performed on account of GO (http://www.geneontology.org). The signi cance of GO term enrichment is indicated standardly by fold changes and P value (FC ≥ 2, P < 0.05 was considered statistically signi cant). 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 that we were interested in. The mRNA related to vascular endothelial proliferation, migration and angiogenesis from the detected mRNA which changed more than two folds with the P value less than 0.05 was chosen. We via cis-acting calculated out the circRNA correlated with the angiogenesis biological process.
GO analysis, which was performed to explore the potential functions of parental gene in terms of the circRNAs, covering three different aspects namely 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)) identi ed by Database for Annotation, Visualization and Integrated Discovery (DAVID; http://www.david.abcc.ncifcrf.gov/).

qRT-PCR validation
qRT-PCR was used to validate the circRNAs expression pattern which were acquired from the microarray detection. Brie y, using the ReverTra Ace qPCR Kit (TOYOBO, FSQ-101) High Capacity cDNA Reverse Transcription Kit with RNase Inhibitor (ABI, 4374966), we reversely transcribed total RNA. qRT-PCR was conducted in 7900 HT Sequence Detection System (ABI, USA) Quant Studio 5 Real-Time PCR System (ABI, USA) using Power SYBR Green PCR Master Mix (ABI, 4368708). We use 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. The circRNAs expression pro les microarray chip assay and data bioinformatics analysis were performed by Shanghai Biotechnology Corporation, Shanghai, People's Republic of China. Based on the expression signal we built a circRNA-miRNA-mRNA interaction network by means of regression model and seed sequence matching analysis. The circRNA-mRNA interaction was predicted with the principle of high related coe cient and at least a perfect seedmatching sequence.

Western blotting analysis
Three days after SCI, followed by rapid isolation of the spinal cords at the lesion site and ash-freezing at -80℃, six mice ( SCI group, n = 3; sham group, n = 3)were euthanized via an overdose of 1% pentobarbital sodium,.The spinal cords were triturated by Tissue Prep (Gering Scienti c Instruments, Beijing, China) and after the designated treatment, cells were dissolved in WB and IP lysis buffer containing phenylmethanesulfonyl uoride (PMSF) and phosphatase inhibitor and centrifuged at 1.3 × 104 g for 15 min at 4 ℃. Subsequently, the protein concentration (3 mg/ml) was quanti ed using a BCA kit (Beyotime Institute of Biotechnology, Jiangsu, China). All samples (20 µl) were separated by SurePAGE (Cat.M00653, GenScript, USA) and 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), rabbit anti-VEGFR2 (Cat# 4108S,ABclonal, 1:1000) were used for probing interested proteins, then transferred onto nitrocellulose membranes that were incubated in primary antibodies at 4 ℃ 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 compared with the β-Actin internal standard.

Statistical analyses
All quantitative data were presented as mean ± standard deviation (mean ± SD). Student's t-test was used to calculate the statistical signi cance. All statistical analyses were performed using SPSS 20.0 software (SPSS, Chicago, IL). The association of the two variables was established using the Pearson product-moment correlation coe cient. The genes that changed more than two folds (including two folds) and P < 0.05 were considered statistically signi cant.

CircRNAs expression pro les in the spinal cord after SCI
The microarrays used in the study contained probes for 37852 mouse circRNAs. Among them, 9364 circRNAs can be found in the circBase and 378 circRNAs were from Deepbase, while the remaining coming from Shbio. As shown in box plot, distributions of the circRNAs expression pro les in all samples were not different (Fig. 1A), indicating that the post-injury changes observed in the levels of individual circRNAs are not random. Differentially expressed circRNAs were sorted by the fold change and p-value (> 2-fold change, p < 0.05), which are shown as red dots and blue dots on the scatter plot ( Fig. 1B) and volcano plot (Fig. 1C). Hierarchical clustering was performed to show the distinguishable circRNAs expression pattern between the two groups (Fig. 1D). The altered circRNAs were clearly separated into two clusters, indicating that the samples had good intra-group consistency, while circRNAs in the SCI group were expressed signi cantly differently from those in the sham group. Compared with the sham group, 1288 circRNAs were observed to be altered (> 2-fold change, p < 0.05) in the spinal cord after SCI, including 991 upregulated and 297 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 and the downregulated circRNAs were transcribed from all chromosomes except chrY ( Fig. 2A).

Real-Time PCR Con rmed circRNAs Expression alteration
In order to verify the microarray data, two circRNAs were selected randomly for real-time PCR validation. Primers used for qRT-PCR of CircRNAs was listed in Table 2. Compared with sham group, the circRNA.7079 were observed to be upregulated (Fig. 2B), while circRNA.12884 were observed to be downregulated signi cantly after SCI (Fig. 2C).

GO analysis and KEGG pathway analysis
We conducted GO analysis and KEGG pathway analysis for the differentially expressed circRNAs according to their host genes to investigate potential biological functions. The GO analysis consequences indicated that biological regulation, cellular process, single-organism process response were the main related biological processes involved in SCI (Fig. 3A). In terms of cellular components, cell, cell part and organelle were highly activated. Most of the molecular function activated in SCI were binding, catalytic activity and molecular function regulator. KEGG analysis showed the top 30 pathways affected by the altered circRNAs that might be involved in the progression of SCI (Fig. 3B). The hostgenes of circRNAs were mainly enriched in the glycosaminoglycan biosynthesis and extracellular matrix(ECM)receptor interaction pathway. The most signi cant pathway was ECM-receptor interaction, which also had the higher enrichment score, re ecting the pathway affected in spinal cord after SCI more than in other diseases or conditions.

CircRNA-miRNA-mRNA ceRNA network prediction and annotation
We screened all the top 20 differently expressed circRNAs after SCI and gured out three of them, circRNA.7079, circRNA.7078 and circRNA.6777, to which related mRNAs could be predicted via ceRNA analysis. The integrated circRNA-miRNA-mRNA network was shown in Fig. 4.

Vascular endothelial proliferation,migration and angiogenesis related circRNAs prediction
We were interested in the circRNAs expression pro les which correlating to vascular endothelial proliferation, migration and angiogenesis after SCI, thus the datasets concerning other functions were excluded. Firstly, to con rm the angiogenesis process related protein expression in spinal cord, western blotting was performed to detect the angiogenesis related protein CD31, VEGF, VEGFR2 expression at 3 days after SCI (Fig. 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, 0.261 ± 0.034 fold of sham, P < 0.001) (Fig. 5B, C, D), revealing the reduction of angiogenesis response. Then, through KEGG analysis, we analyzed the pathways correlating to endothelial proliferation, migration and angiogenesis and counted the numbers of related circRNAs (Fig. 5E). The host-genes of circRNAs were mainly enriched in the ECMreceptor interaction (number = 97) and PI3K-Akt signaling pathway(number = 85) as shown in Fig. 5E. As for in VEGF signaling pathway, the host-genes of circRNAs, such as PLCy, SPK, Ras, COX2, MAPKAFK, HSP27 were upregulated, while PKC, NEAT were downregulated (KEGG mmu04370, Figs. 5F). VEGFR2 is the main receptor of VEGFA to trigger a series of vascular endothelial proliferation, migration and angiogenesis through VEGF signaling pathway as shown in the KEGG mmu04370.

Discussion
Studies have described that altered levels of speci c circRNAs are involved in many central nervous system diseases, such as brain injury, ischemia, stroke, cancer, neurodegenerative diseases [14][15][16][17][18][19]. So far, the functions of most of the SCI altered circRNAs have not been determined. In this study, we have comprehensively elucidated the expression pro les of circRNAs in a mouse model of SCI.
We screened all the top 20 differently expressed circRNAs and gured out three of them, circRNA.7079, circRNA.7078 and circRNA.6777, to which related mRNAs could be predicted via ceRNA analysis. The fold change of circRNA.7079 was the TOP 1, and the alteration was con rmed by qRT-PCR. It was interesting to nd that circRNA.7079 was involved in the regulation of mRNA of BIRC5, Lgals3 and Pbk. Moreover, ceRNA analysis of circRNA-miRNA-mRNA predicted that circRNA.7079 may affect the expression of target mRNAs by sponging mmu-miR-326-5p, mmu-miR-761, mmu-miR-1907 or mmu-miR-3473. It was notable that one of the target mRNA,Birc5 can promote the expression of VEGF and regulate the expression of angiogenesis-associated factor and increase the angiogenesis and migration [20]. In addition, extensive studies of Lgals3 have shown that it is expressed in several tissues acting as a regulator of cell proliferation, apoptosis, adhesion, and migration, then affecting many aspects of events, such as angiogenesis and tumorigenesis. [21]. Besides, circRNA.7078 may regulate the expression level of Top2a, Cd44, Ms4a6c and Knstrn by sponging mmu-miR-761, mmu-miR-1907 or mmu-miR-3473 and mmu-miR-5710. It was reported that Top2a was the top 10 core genes in a SCI model at 2 days after spinal cord injury [22]. Top2a is expressed in proliferating cells during the S and G2M phases of the cell cycle and is primarily involved in regulating cell proliferation [23].Moreover, circRNA.6777 may regulate the expression level of Fn1 by sponging mmu-miR-7684-5p. Extensive studies of Fn1 have shown that it inhibits apoptosis and regulating EMT to promote melanoma proliferation and metastasis [24]. Recently, it was indicated that down-regulation of FN1 suppresses proliferation, migration and invasion to inhibit colorectal carcinogenesis [25].
Angiogenesis is central to a number of pathophysiological process, from embryogenesis to wound healing [26]. It is a complex process including extensive interplay between endothelial cells (EC), soluble factors and extracellular matrix components. It was reported that vascular-related events such as vascular endothelial proliferation, migration and angiogenesis promote the recovery of SCI [27]. CD31(PECAM-1) involved in angiogenesis and the interactions of endothelial cell-cell adhesion molecules are important in the formation of new vessels, according to the previous studies [28]. As a main angiogenetic factor, vascular endothelial growth factor A (VEGFA) promotes endothelial cell proliferation after injury in various types of tissue, including the central nervous system (CNS) [29]. It is reported that VEGF signaling pathway initiated by VEGFA/VEGFR2 leads to endothelial cell proliferation, migration, survival and new vessel formation involved in angiogenesis [30]. It is worth taking note that endothelial cells start newly forming microvessels from 3 days after injury, then restoring the density of microvessels to a normal level by one week after SCI [31]. In our present studies, western blot analysis con rmed a decreased protein expression of CD31, VEGF and VEGFR2 at 3 days after SCI. Therefore, we supposed that microvessels were damaged during the early 3 days after SCI, but at the same time the endogenous signaling pathways were gradually activated to initiate vascular endothelial proliferation, migration and angiogenesis. Then, based on KEGG pathways analysis, we found there are 10 endogenous signaling pathways associated with vascular endothelial proliferation, migration and angiogenesis. The most altered circRNAs were enriched in the ECM-receptor interaction and PI3K-Akt signaling pathway. We are especially interested in the host genes of altered circRNAs in VEGF signaling pathways associated with vascular endothelial proliferation, migration and angiogenesis after SCI. We found that there are 8 host genes of circRNAs such as PLCy, PKC, SPK, Ras, NFAT, COX2, MAPKAFK, HSP27 were involved in the VEGF signaling pathway.
So far, some circRNAs that have been reported to play role in vascular endothelial proliferation, migration and angiogenesis. For example, circRNA hsa_circ_0003575 induces vascular endothelial cells proliferation and angiogenesis by regulating oxLDL [32]. It is reported that up-regulating circRNA-MYLK promoted the growth, angiogenesis and metastasis of BC xenografts [33]. A recent study indicated that circN x downregulation promoted CM proliferation and angiogenesis and inhibited CM apoptosis after MI, attenuating cardiac dysfunction and improving the prognosis [34]. However, the circRNAs involved in vascular endothelial proliferation, migration and angiogenesis after SCI are not being clari ed. We attempted to correlate mRNAs and circRNAs by ceRNA prediction in order to determine the vascular endothelial proliferation, migration and angiogenesis related circRNAs based on the GO and KEGG results. We found out 93 circRNAs were involved in vascular endothelial proliferation, migration and angiogenesis. In summary, these ndings may be helpful to understand the mechanisms of the early stage of SCI, and might present novel molecular targets for clinical therapy of SCI.

Conclusion
This abundant datasets suggest that circRNAs locally regulate their related protein-gene expression and play key roles in the vascular endothelial proliferation, migration and angiogenesis of SCI. Our ndings provide a potential new possibility for the treatment of SCI in mice by modulating circRNAs.

Declarations
Ethics approval and consent to participate The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. All experimental procedures were carried out according to the World Medical Association Declaration of Helsinki and by the guidelines and regulations approved by the Zhejiang University School of Medicine.

Consent for publication
Written informed consent for publication was obtained from all participants.

Availability of data and materials
The datasets used or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests
All authors con rm that there is no competing interests regarding this study.                    factor on the X-axis indicates the speci city of a pathway activating in SCI samples compared with sham samples.
The bubble size represents the number of circRNAs enriched in a pathway. Red bubble indicates more signi cant pathway.

Figure 4
The biomathematical about ceRNA analysis of circRNA-miRNA-mRNA crosstalk network in spinal cord of mice after SCI. circRNA.7079, circRNA.7078 and circRNA.6777, were plotted as green rhombi. And miRNAs targeted by circRNAs on the network were plotted as red arrow. mRNAs that might be regulated by these three circRNAs were plotted as yellow dots. The biomathematical about ceRNA analysis of circRNA-miRNA-mRNA crosstalk network in spinal cord of mice after SCI. circRNA.7079, circRNA.7078 and circRNA.6777, were plotted as green rhombi. And miRNAs targeted by circRNAs on the network were plotted as red arrow. mRNAs that might be regulated by these three circRNAs were plotted as yellow dots.

Figure 4
The biomathematical about ceRNA analysis of circRNA-miRNA-mRNA crosstalk network in spinal cord of mice after SCI. circRNA.7079, circRNA.7078 and circRNA.6777, were plotted as green rhombi. And miRNAs targeted by circRNAs on the network were plotted as red arrow. mRNAs that might be regulated by these three circRNAs were plotted as yellow dots.
Page 30/49 Figure 4 The biomathematical about ceRNA analysis of circRNA-miRNA-mRNA crosstalk network in spinal cord of mice after SCI. circRNA.7079, circRNA.7078 and circRNA.6777, were plotted as green rhombi. And miRNAs targeted by circRNAs on the network were plotted as red arrow. mRNAs that might be regulated by these three circRNAs were plotted as yellow dots.