The successful establishment of rat SCI model, the isolation of blood exosomes, and the identification of exosomes were detailedly described in our recently published paper (Li et al., 2023). In brief, six adult female rats were randomly divided into two groups, the experimental group (EG) in which SCI surgery was performed on rats, and the control group (CG) in which the skin was injured on rats. Both cohorts of rats were maintained under the same feeding conditions before and after surgery. After 24 h, blood was collected and exosomes were extracted by using the density gradient centrifugation. Then, the classical characteristics of isolated exosomes were verified through Western blot, electron microscopy, and nanoparticle tracking analysis (NTA). Related data on this part would be referred to our previous publication (Li et al., 2023).

According to the manufacturer’s protocol, Trizol reagent (Thermo Fisher Scientific, Waltham, MA, USA) was used to obtain high-quality RNA from circulating exosomes. Concentrations of different RNA samples were measured through using Nanodrop Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Then the quantitative control as well as the integrity control were performed for isolated RNA prior to experiments in the next step.

In the process of quantitative control, the purity of RNA was evaluated by using OD260/OD280 values. The final results would be considered as “Pass”, if these OD260/OD280 values are within the range between 1.8 and 2.1, suggesting isolated RNA has high quality and can be used further.

Regarding the integrity control, 1% agarose gel (Thermo Fisher Scientific, Waltham, MA, USA) was utilized to assess the integrity of isolated RNA in this study.

Prior to RNA library construction, ribosomal RNA (rRNA) needs to be completely removed from total RNA. To this end, NEBNext rRNA Depletion Kit (New England Biolabs, Ipswich, MA, USA) was used according to the producer’s protocol. Following this, the RNA library was then constructed through using NEBNext^® Ultra™ II Directional RNA Library Prep Kit (New England Biolabs, Ipswich, MA, USA) according to the standard procedures. Furthermore, RNA library was quality-controlled and then quantified via a BioAnalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). In the end, circRNA sequencing was conducted on an illumina HiSeq sequencer from Cloud-Seq Biotech (Shanghai, China). All raw data were uploaded to GEO repository as the GSE213561 study (GSM6589180; GSM6589181; GSM6589182; GSM6589183; GSM6589184; GSM6589185), and the link is as follows: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE213561.

As described in the above part, raw data which contain original paired-end reads were initially acquired from Illumina HiSeq 4000 systems (Cloud-Seq Biotech, Shanghai, China). Through using Cutadapt software,¹ 3′ adaptor was modified, and low-quality reads were removed in order to keep high-quality reads for further analysis (Han et al., 2020). To display relative levels of these circulating circRNAs, these modified reads were analyzed. HISAT2 software² was used to make an alignment of these pre-selected reads to the human reference genome. Next, according to gene transfer format (GTF) files which were obtained from the Ensembl database,³ the expression levels of circRNAs, characterized by the fragments per kilobase of exon model per million mapped fragments (FPKM), were calculated via Cuffdiff software (Trapnell et al., 2010). Following this, specific fold changes and corresponding P-values were calculated based on the FPKM index, and differentially expressed circRNAs from circulating exosomes were eventually determined.

Moreover, in order to predict their probably involved pathways and functional processes, KEGG pathway analysis⁴ and GO analysis⁵ were performed for circulating upregulated as well as downregulated circRNAs, respectively (Kanehisa and Sato, 2020). The value –log10(P-value) was calculated and presented as the enrichment score. For both KEGG pathway and GO analysis, fold change > 2 and P-value < 0.05 were taken as the threshold of the differential expression of these circRNAs.

Next, we would like to separately verify the top ten dysregulated circulating circRNAs in SCI rats, which were revealed by RNA sequencing previously. To this end, RT-PCR was utilized here to show the relative levels of these 10 circulating exosomal circRNAs from both group rats. As mentioned in the above part, Trizol reagent was used to isolate RNA, and then cDNA was synthesized by utilizing SuperScript™ IV First-Strand synthesis kit (Thermo Fisher Scientific, Waltham, MA, USA). For RT-PCR process, SYBR Green master mix (CloudSeq, Shanghai, China) was additionally employed. Primer 5.0 software (PREMIER Biosoft, Palo Alto, CA, USA) was used to design all primer sequences for circRNA targets as well as the housekeeping gene, i.e., glyceraldehyde-3-phosphate dehydrogenase (GAPDH). All primers were purchased from CloudSeq Biotechnology (Shanghai, China). The primer sequences are shown in Table 1. All raw data were acquired by QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA), and further 2^–ΔΔCt method was applied to show the relative expression of these circulating circRNAs. Given that most of differentially expressed circRNA are novel, we used their chromosome locations to name them.

Considering the importance of functional interactions among circRNAs, mRNAs and miRNAs during the gene regulation, the upstream as well as the downstream on the genome of these differential expressed circRNAs were extended in order to find functional genes near these circRNAs. To this end, corresponding miRNAs were predicted through employing miRNA target prediction software miRanda⁶ based on these differentially expressed circRNAs, and then potential mRNAs were further anticipated by applying Targetscan⁷ as well as miRDB⁸ (Sun et al., 2020). In the end, circRNA-miRNA-mRNA interaction network was established and visualized via Cytoscape⁹ (Zhao et al., 2021), where these differentially expressed circRNAs from circulating exosomes were centered.

GraphPad Prism 8 software was utilized to make statistical analyses for these experimental results. Data were shown as means ± standard deviation (SD) from more than three statistical independent experiments. To compare the statistical significances among different groups, t-test of two independent samples was utilized, respectively. When P-value < 0.05, the difference was regarded as statistically significant. Especially, *P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001.

The experimental design of this study was shown in Figure 1. Part data of preliminary work including SCI rat model establishment and exosome identification were presented in our recently published paper (Li et al., 2023).

Through enrichment analyses, a total of 709 types of circRNAs were found to be downregulated whereas 346 kinds of circRNAs were shown to be upregulated in SCI rats in comparison with the control rats. In this experiment, three SCI rats and plus three control rats were included. Through the normalization as FPKM, the heat map of differentially expressed circRNAs from circulating exosomes including 15 kinds of downregulated circRNAs and 7 kinds of upregulated circRNAs was especially visualized (Figure 2A), suggesting that there exists a significantly different expression pattern of circulating exosomal circRNAs between SCI rats and the control rats (P < 0.05). In a similar vein, the volcano plots demonstrated consistent changes of differentially expressed circRNAs in SCI rats (Figure 2B). Taken together, these data give the first impression that some circulating circRNAs display differential expression in SCI rats.

To further speculate their specific functional activities and involved signaling pathways in the context of SCI, both GO analysis and KEGG analysis were next carried out in this study. Figures 3, 4 together demonstrated the top enriched GO terms and KEGG pathways for these differentially expressed cicRNAs. Especially, upregulated circRNAs were analyzed in Figure 3, and downregulated circRNAs were analyzed in Figure 4. These results together provide an overview of these probably involved cellular components, biological processes, and molecular functions as well as signaling pathways for differentially expressed circRNAs in SCI rats (Figures 3, 4), which may offer some meaningful indications for future research.

The top five enriched biochemical pathways for upregulated circRNAs were African trypanosomiasis, Malaria, non-small lung cancer, small cell lung cancer, and purine metabolism, as demonstrated by KEGG pathway analysis based on located genes (Figure 3A). On the other hand, GO enrichment analysis showed that these upregulated circRNAs were enriched in plenty of biological processes associated with inflammation, cellular components, and molecular functions. The top five enriched terms of biological process were oxygen transport, negative regulation of proteolysis involved in cellular protein catabolic process, negative regulation of proteolysis, negative regulation of protein catabolic process, and negative regulation of proteasomal ubiquitin-dependent protein catabolic process (Figure 3B). The top four enriched terms of cellular component were hemoglobin complex, cytosolic part, cytosol, and cytoplasmic part (Figure 3C). At last, oxygen transporter activity, oxygen binding, heme binding, tetrapyrrole binding, and transition metal ion binding were the top five enriched terms of molecular function (Figure 3D). Of special interest, the potential involvement of some circRNAs in tumors may indicate the potency of circRNA-elicited microenvironment changes in inflammatory processes.

The top eight enriched biochemical pathways were GABAergic synapse, ErbB signaling pathway, mRNA surveillance pathway, T cell receptor signaling pathway, glutamatergic synapse, thyroid hormone signaling pathway, axon guidance and RNA transport, as demonstrated by KEGG pathway analysis based on located genes (Figure 4A). In addition, GO enrichment analysis showed that these downregulated circRNAs were enriched in plenty of biological processes associated with inflammation, cellular components, and molecular functions. The top five enriched terms of biological process were transcription from RNA polymerase II promoter, RNA metabolic process, regulation of transcription, DNA-templated, and positive regulation of T cell activation (Figure 4B). The top five enriched terms of cellular component were SWI/SNF superfamily-type complex, SWI/SNF complex, proteasome complex, pericentric heterochromatin, and npBAF complex (Figure 4C). At last, transcription factor binding transcription factor activity, transcription factor activity, RNA polymerase II activating transcription factor binding, protein binding transcription activity, and protein-lysine N-methyltransferase activity were the top five enriched terms of molecular function (Figure 4D). In combination with data for upregulated circRNAs, these bioinformatics data together suggest there are different pathways and functions involved for these dysregulated circRNAs.

As addressed in the above part, several significantly upregulated and downregulated circRNAs were preliminarily selected through RNA sequencing (Figure 2). Later on, the potential roles of these circulating circRNAs in SCI rats were predicted via KEGG and GO analyses (Figures 3, 4). Therefore, we would like to take advantage of qRT-PCR assay, to verify whether these candidate circRNAs are significantly dysregulated indeed. It turned out that there were a total of ten kinds of differentially expressed circRNAs, which were further verified by qRT-PCR in this study. Quantification results displayed that the relative expression of chr4:208359914–208362182+, chr15:20088296–20092102+, chr1:175098934– 175134845–, chr1:175099657– 175128203–, and chr1:175104454– 175134845–, was significantly elevated in SCI rats in comparison with the control rats (Figure 5A). By contrast, the relative levels of chr11:74154652– 74159524–, chr12:45412398– 45412635–, chr7:137630261– 137648924–, chr6:6280974–6281188+, and chr4:225251864–225254087+ were dramatically diminished in SCI rats compared to the controls (Figure 5B). Moreover, the detailed information of these ten candidate circRNAs is summarized in Table 2. Given that most of differentially expressed circRNA are novel and the gene symbol of three kinds of circRNA is the same, we named them according to their chromosome locations. Relevant predicted pathway information of each circRNA can be checked through circBase¹⁰ (Table 2). These quantitative results not only provide some promising targets for future in vivo and in vitro research, but also validate the expression profile of circulating circRNAs indicated by RNA sequencing.

Based on ceRNA mechanism, we would like to take advantage of circRNA-miRNA-mRNA interaction network to anticipate potential functions of circRNAs in the injured spinal cord. As shown in Figure 6, five kinds of upregulated circRNAs as well as five kinds of downregulated circRNAs verified by RT-PCR were chosen and centered in this network. These circRNAs could co-express and interact with 37 miRNAs which are indicated by red triangles and 145 mRNAs which are indicated by blue boxes (Figure 6). All these upregulated circRNAs had four to six target miRNAs, as exemplified by chr4:208359914–208362182+ which has five target miRNAs including rno-miR-26b-3p, rno-miR-17-1-3p, rno-miR-3068-5p, rno-miR-3564, and rno-miR-20b-3p (Figure 6). In the meanwhile, these downregulated circRNAs such as chr11:74154652– 74159524–, chr12:45412398– 45412635–, chr7:137630261– 137648924–, chr6:6280974–6281188+, and chr4:225251864–225254087+ had three to six target miRNAs generally (Figure 6). Thus, these bioinformatics data together imply the potential important functions of the above-mentioned circulating circRNAs in regulating other genes.