Role of circular RNA expression in the pathological progression after spinal cord injury

Wen-Zhao Wang, Jun Li, Lei Liu, Zheng-Dong Zhang, Ming-Xin Li, Qin Li, Hui-Xu Ma,Hai Yang, Xiao-Ling Hou

Abstract Differential expression of non-coding RNA after traumatic spinal cord injury (TSCI) is closely related to the pathophysiological process. The purposes of this study were to systematically profile and characterize expression of circular RNA (circRNA) in the lesion epicenter of spinal tissues after TSCI, and predict the structure and potential function of the regulatory circRNA/miRNA network. Forty-eight C57BL/6 mice were randomly and equally assigned to two groups: one subjected to TSCI at T8-10 with an Allen’s drop impactor, and a second subjected to laminectomy without TSCI. Spinal cord samples were stained with hematoxylin and eosin, sequenced, and validated. RNA-Seq, Gene Ontology analysis, Kyoto Encyclopedia of Genes and Genomes analysis, and network analyses (Targetscan and miRanda) were used to predict and annotate the circRNA/miRNA/mRNA network. Luciferase reporter, quantitative reverse transcription polymerase chain reaction, and western blot assays were used to profile expression and regulation patterns of the network in mouse models of TSCI. Hematoxylin-eosin staining revealed severe damage to the blood-spinal cord barrier after TSCI. Differentially expressed circRNA and miRNA profiles were obtained after TSCI; differentially expressed circRNAs, which were abundant in the cytoplasm, were involved in positive regulation of transcription and protein phosphorylation. miR-135b-5p was the most significantly downregulated miRNA after TSCI; circRNAAbca1 and KLF4 were predicted to be its target circRNA and mRNA, respectively. Subsequently, the circAbca1/miR-135b-5P/KLF4 regulatory axis was predicted and constructed,and its targeted binding was verified. After inhibiting circAbca1, GAP43 expression was upregulated. Differential expression of circRNAs might play an important role after TSCI. circAbca1 plays a neuroinhibitory role by targeted binding of the miR-135b-5P/KLF4 axis. The identified circRNA/miRNA/mRNA network could provide the basis for understanding pathophysiological mechanisms underlying TSCI, as well as guide the formulation of related therapeutic strategies. All animal protocols were approved by the Research Ethics Committee of West China Hospital of China (approval No. 2017128) on May 16, 2017.

Key Words: bioinformatics; circRNA/miRNA/mRNA network; circular RNA; gene; KLF4; miR-135b-5p; spinal cord injury; trauma Chinese Library Classification No. R452; R363; R364

Introduction

Traumatic spinal cord injury (TSCI) is a type of central nervous system injury with high mortality and morbidity (Chen et al.,2020; Jiang et al., 2020; Ray et al., 2020). Damage often leads to loss and death of nerve tissues, causing sensory and motor dysfunction (Fouad et al., 2011; Holmes, 2017). Motor vehicle collision is the most common cause of TSCI. However, with the global population aging, the proportion of TSCI due to falls has increased from 16% in 2000 to 30% in 2012 (Cristante et al.,2012). Peripheral nerve fibers project from the spinal cord to relay signals to and from muscles and organs. However, unlike peripheral nerves, neural connections in the spinal cord are rarely repaired after injury; accordingly, almost all patients with TSCI experience pain throughout their lifetime. Current understanding of the underlying mechanisms of TSCI remain limited and there is no effective treatment for this condition(Tran et al., 2018). Therefore, it is necessary to better understand the cellular and molecular pathophysiological mechanisms following TSCI to identify new therapeutic strategies that enhance axonal regeneration and functional recovery.

Modern biological techniques have confirmed that 98%of the human genome comprises RNAs that do not code for protein - most are called noncoding ribonucleic acids(ncRNAs). ncRNAs play key roles in the regulation of genetic programming, epigenetics, transcription, and translation in many pathophysiological processes (Geng et al., 2018).Circular RNAs (circRNAs), a recently discovered subtype of the ncRNA family, have a closed-loop structure with no 5′to 3′ polarity or polyadenylated tail. Compared with linear RNAs, circRNAs have a more stable structure and are highly expressed (Zhang and Xin, 2018). Studies have found that circRNAs such as circRim2, circEl2, circPlxnd1, and circDym are abundant and function in the nervous system (Rybak-Wolf et al., 2015; Wang et al., 2018a). However, further investigation is required to predict and verify the functions of many newly identified circRNAs. MicroRNAs (miRNAs) are evolutionarily conserved, single-stranded ncRNAs composed of 18-23 nucleotides. miRNAs regulate gene expression at the post-transcriptional level, either by degradation or inhibiting translation, by binding to the 3′ UTR of target mRNAs (Ameres and Zamore, 2013). circRNA/miRNA/mRNA networks have been shown to play important regulatory roles in several studies (van Rossum et al., 2016; Han et al., 2018; Kleaveland et al., 2018). However, roles of the circRNA/miRNA/mRNA network in TSCI remain to be systematically investigated.Our present study employed high-throughput sequencing of miRNAs and circRNAs to screen for significantly expressed miRNAs following TSCI. After matching these miRNAs with differentially expressed circRNA profiles, signaling pathway and function analysis were performed, and the resulting predictions were used to construct an interaction network.After verifying the predicted targeted binding relationships,we preliminarily validated the interaction axes in vitro.

The aim of this study was to predict and construct a circRNA/miRNA/mRNA network that plays a potential regulatory role after TSCI.

Materials and Methods

Animals

All animal protocols were approved by the Research Ethics Committee of West China Hospital, China (No. 2017128) on May 16, 2017. Animal experiments met the requirements of the ARRIVE Guidelines (https://arriveguidelines.org/). Fortyeight adult male mice provided by the Experimental Animal Center of Sichuan University (license No. SCXK (Chuan) 2018-026), specific pathogen-free, 8 weeks of age, body weight 18-22 g, were used in animal experiments. Animals were housed in a clean facility under controlled conditions of 22-24°C,relative humidity of 30-50%, and a 12-hour light/dark cycle.

Construction of a mouse SCI model and tissue acquisition

Mice were randomly divided into six groups (eight mice per group): three underwent SCI (SCI1, SCI2, and SCI3),and the remaining three underwent a sham operation with laminectomy only (Sham1, Sham2, and Sham3). Mouse models were established as previously described (Wang et al.,2018c). Briefly, after anesthesia with pentobarbital sodium(30 mg/kg), tissue layers were incised, exposed, and localized to the thoracic T8-10 segment. In sham groups, vertebral plates were removed before suture. In SCI groups, contusion injuries were created using an Allen’s drop impactor with a drop weight of 6 g and drop height of 60 mm. At the moment of injury, the mouse tail was lifted high or swung violently, and the lower limbs convulsed violently, which were regarded as criteria for the successful induction of injury induction. The wound was then closed by suturing (Xie et al., 2018). Tissues of the lesion epicenter were extracted from seven randomly selected mice per group and combined as one sample for sequencing and qRT-PCR. Tissues from the lesion epicenter of the remaining mice in each group were used for hematoxylin and eosin (HE) staining. Spinal cord samples were collected separately after 5 days. After incision, spinal cord tissue in the lesion epicenter was dissected and subjected to RNA extraction by adding TRIzol (Invitrogen, Carlsbad, CA, USA),followed by chloroform. The middle layer was extracted after centrifugation, and isopropanol was added. The resulting mixture was centrifuged and the supernatant was discarded.Ethanol was added and the mixture was further centrifuged,followed by drying of the precipitate. RNase-free water was added to dissolve the resulting pellet.

HE staining of tissues at the site of spinal cord injury

At 5 days post-SCI, mice were administered anesthesia and a previously described protocol was used to harvest and stain spinal cord tissues (Wang et al., 2018c). Briefly, 4%paraformaldehyde was perfused into mice for fixation, and the lesion epicenter of spinal tissues was extracted and placed in 4% paraformaldehyde overnight, followed by dehydration,paraffin embedding, and sectioning. After xylol dewaxing,samples were dehydrated in an ethanol gradient, stained with hematoxylin, differentiated with 1% hydrochloric acid,and stained with eosin. After dehydration using an ethanol gradient, samples were immersed in xylol and mounted in neutral gum. Tissue integrity was compared between SCI and sham groups under a light microscope (Olympus Corporation,Tokyo, Japan).

Sequencing library construction and sequencing analysis

Total RNA was extracted from spinal cord tissue and measured with a spectrophotometer (Invitrogen) according to the ratio of optical densities at 260 nm and 280 nm (OD260/OD280). To generate a sequencing library for circRNA, total ribosomal RNA was first removed. Purified and recovered residual RNA was randomly broken into short fragments using fragmentation reagents. Fragmented circRNA was used as a template to synthesize single-stranded cDNA with a six-base random primer, before further synthesis into double-stranded cDNA using dNTPs, RNaseH, and DNA polymerase I. Double-stranded products were purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA), and the sticky ends of DNA were repaired to blunt ends using T4 and Klenow DNA polymerases.After adding adenosine at the 3′ end, a linker was added.Subsequently, AMPureXP beads were used for fragment selection. Finally, PCR amplification was performed to establish a sequencing library. The Hiseq2000/2500 platform(Illumina, San Diego, CA, USA) was used for sequencing with a read length parameter of 2 × 150 bp. The miRNA library was treated with a TruSeq Small RNA Sample Prep Kit (Illumina)and sequenced with an Illumina Hiseq 2000/2500. Library sequencing was performed as single-ended 1 × 50 bp reads.

Expression profile analysis and differential expression data collation

After obtaining the original sequencing data, differential expression was identified and analyzed. RNA circularization can be divided into exon circularization and intron circularization, both have reverse cleavage sites that are used to identify the sequence of the circRNA. CIRCExplorer was used to predict the circRNA sequence, and TopHat was used for alignment. circRNA positional relationship classification was used for annotation, and the distribution of gene expression values was calculated and compared between different samples to confirm their biological reproducibility.Finally, statistical analysis of differential gene expression in different tissue samples was performed. The coincidence of circRNA expression in SCI and Sham groups was analyzed using a Venn diagram, and circRNAs with low expression were excluded using the following conditions:P＜ 0.05 and|log2(fold change)| ≥ 1. At-test was used for comparison, and differential expression factors were arranged in descending order. Selected miRNA sequences were compared with those of known miRNAs in the miRBase 21.0 database (http://www.mirbase.org/) and used for RNA secondary structure prediction. Finally, miRNAs of different tissue samples were statistically compared, and those that were significantly expressed after injury were selected.

Gene Ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of differentially expressed circRNA

The results of differential expression analysis were further analyzed for function and signaling pathway enrichment. In GO enrichment analysis (Gene Ontology, 2015), the basic unit is the term, with each term corresponding to an attribute.GO enrichment analysis first maps all significant differentially expressed genes to terms in the GO database (http://www.geneontology.org/), and then calculates the number of genes for each term and employs a hypergeometric test. Enriched genes were compared with the whole genome background withP≤ 0.05 as the threshold. GO terms that were significantly enriched among the significantly differentially expressed genes were screened, and the main biological functions of differentially expressed genes were predicted.GO term enrichment analysis results represent biological processes, cellular components, and molecular functions.KEGG (http://www.genome.jp/kegg/) (Kanehisa et al., 2017)the main public database of pathways, is commonly used for pathway enrichment analyses. This hypergeometric test identifies significantly enriched pathways, representing the most important biochemical metabolic and signal transduction pathways associated with differentially expressed genes.

miRNA/circRNA interaction prediction analysis

The results of differential expression analysis were comprehensively analyzed to predict targeted interaction relationships. Targetscan (http://www.targetscan.org/) and miRanda (http://www.microrna.org/) (Zhang et al., 2014)software were used to predict miRNA target genes and miRNA/circRNA interactions. Targetscan predicts complete sequence complementation information by seed matching in the seed region. miRanda constructs a scoring matrix to predict target information by base complementation (A=U,G=C). Default parameters were used in Targetscan, and a Max_Energy threshold ＜ -10 was set during miRanda analysis.Overlapping results obtained by these two algorithms were considered the final result of this interaction analysis.Differential circRNA analysis was performed using the default difference threshold, |log2(fold change)| ≥ 1, andP≤ 0.05.According to the number of reverse-cleaved back-splicing sites of the circRNA, accuracy was high and the proportion of false positives was low. The threshold for the TargetScan value was set to 50, such that the higher the value, the higher the combination rate. The threshold set was -10 for miRanda, as smaller values allow for higher combination rates. Finally, a miRNA/circRNA interaction map was drawn using Cytoscape to visualize the predicted interaction network information.

qRT-PCR quantitation

To verify sequencing results, total RNA extracted fromtissue was subjected to PCR detection. qRT-PCR was used to quantitate expression ofcircRNAAbca1(circAbca1),mmumiR-21a-5p_R+1,KLF4,GAPDH, andU6(Table 1). A reverse transcription reagent kit (TUREscript 1st Stand cDNA SynthesisKit, Aidlab, Beijing, China) was used for reverse transcription,and cDNA was used as a template to prepare a 20-µL reaction system for PCR in an ABI PRISM 7500 RT-PCR System (Applied Biosystems, Foster City, CA, USA). SYBR green dye was used to detect amplification in real time. GAPDH was used as an internal reference for circRNA and mRNA, while U6 was used as an internal reference for miRNA expression analysis. Predenaturation at 95°C was performed for 3 minutes, followed by 40 cycles of 95°C for 10 seconds and annealing/extension at 58°C for 30 seconds. miRNA amplification conditions were pre-denaturation at 95°C for 10 minutes, followed by 40 cycles of 95°C for 2 seconds and annealing/extension at 60°C for 30 seconds. The experiment was repeated three times and data were analyzed using the 2-∆∆Ctmethod (Wang et al.,2018c).

Table 1 ｜Primers used in qRT-PCR

Luciferase reporter assay

To verify the predicted targeted binding relationship, 293T cells were seeded into 96-well plates at 70% confluence and transfected with the experimental plasmid pMIR-REPORT Luciferase-circAbca1(WT) (LC Science, Houston, TX, USA)or pMIR-REPORT Luciferase-circAbca1(MUT1 + MUT2) (LC Science) after 24 hours. Subsequently, a luciferase reporter plasmid (Firefly andRenilla) (Invitrogen) was transfected with 0.25 µL of microRNA transfection reagent (Invitrogen) for 48 hours in each well. A dual-luciferase reporter assay (E1910;Promega, Madison, WI, USA) was then performed with 50 µL of pre-mixed Stop & Glo Reagent (Invitrogen). Results were recorded at 2-second intervals.

Cell culture and transfection

The targeted binding relationship was preliminarily verifiedin vitro. Primary C57BL/6 mouse dorsal root ganglion cells(PriCells, Wuhan, China) were cultured at 37°C and 5% CO2in Neurobasal medium (Thermo Fisher Scientific, Waltham, MA,USA) supplemented with 1% B-27TMserum-free supplement(50X, Thermo Fisher Scientific) and 100 kU/L of penicillin and streptomycin (Thermo Fisher Scientific). Cultured cells were transfected with circAbca1 inhibitor (5′-CCA ACT CAA TGT CAG CTG T-3′) or the corresponding controls (Ribobio, Guangzhou,China). Next, 120 µL of 1× riboFECTTMBuffer was diluted with 5 µL and 20 µM siRNA stock solution, and then 12 µL riboFECTTMReagent (Ribobio, Guangzhou, China) was added and mixed at room temperature for 10 minutes. The mixture was added to medium without penicillin and streptomycin,and then added to cells for an additional 48 hours of culture.

Western blot analysis

After transfection, cell samples were collected and total protein was extracted. Proteins were incubated with cold lysis buffer (Solarbio, Beijing, China), centrifuged,concentrated, electrophoresed, and subsequently transferred onto polyvinylidene fluoride membranes. After blocking with phosphate-buffered saline containing 5% non-fat milk overnight at 4°C (Wang et al., 2018c), membranes were incubated with antibodies against growth-associated protein 43 (GAP43; rabbit monoclonal, 1:1000; Abcam, Cambridge,UK) and GAPDH (rabbit monoclonal, 1:2500; Abcam)overnight at 4°C. The following day, membranes were treated with horseradish peroxidase-conjugated secondary antibodies(goat anti-rabbit IgG, 1:5000; Proteintech, Wuhan, China). An enhanced West Pico ECL Substrate kit (Solarbio) was utilized to perform chemiluminescent detection. GAP43 expression was quantified using ImageJ software (National Institutes of Health, Bethesda, MD, USA) and normalized to GAPDH.

Statistical analysis

Statistical analysis was performed using Prism 8 (GraphPad Software, La Jolla, CA, USA) and SPSS 20.0 software (IBM,Armonk, NY, USA). Data are expressed as mean ± standard deviation (SD). Two-samplet-test, one-way analysis of variance, and Dunnett’s test were used to assess statistical significance.P-values ＜ 0.05 were considered statistically significant.

Results

Establishment of SCI mouse models and expression profiles of circRNAs

To explore expression profiles of circRNAs in the lesion epicenter after SCI, a mouse model of SCI was established using an Allen’s drop method. HE-stained sections of SCI and sham groups were observed by microscopy. Nuclei were stained dark blue by hematoxylin, while cytoplasm appeared pink and collagen fibers were pale pink after eosin staining.The spinal cord structure in sham groups was clear with intact and unruptured meninges, and a centered central tube. Moreover, there was no bleeding or signs of damage(Figure 1A). In contrast, the spinal cord structure of SCI groups exhibited significant damage and rupture, meningeal defects,and central tube offset deformation. Extensive red blood cell diffusion and inflammatory cell infiltration caused by vascular rupture were also observed (Figure 1B). A cluster map and volcano plot show differentially expressed circRNAs after TSCI(Figure 1CandD).

Enriched GO terms and KEGG pathways of differentially expressed circRNAs

To explore major biological functions and significantly enriched pathways, GO and KEGG analyses were performed.Differential expression profiles of circRNAs were further screened according to their expression level and log2(fold change). Enriched terms for differentially expressed genes identified by GO term analysis were used to determine their major biological functions. The results of GO enrichment analysis showed that the highest degree of enrichment was associated with ‘negative regulation of plasminogen activation’, ‘peptide cross linking’, ‘positive regulation of the endothelial cell apoptotic process’, and ‘macropinosome’(Figure 2A). Analysis of significantly enriched KEGG pathways revealed the most important biochemical metabolic pathways and signal transduction pathways encoded by differentially expressed genes included ‘complement and coagulation cascades’, ‘glycerolipid metabolism’, ‘one carbon pool by folate’, ‘pentose phosphate pathway’, ‘bladder cancer’, and‘P53 signaling pathway’ (Figure 2B).

Prediction of interactions between differentially expressed miRNA/circRNAs

To investigate regulatory roles of the targeted interaction axis after SCI, we predicted the miRNA/circRNA axis using bioinformatic analyses. After TSCI, differential expression profiles of circRNAs and miRNAs were screened according to their expression and log2(fold change).miR-135b-5pwas significantly downregulated (P＜ 0.05), as it satisfied the significance threshold of |log2(fold change)| ≥ 1. Indeed,miR-135b-5pshowed the most significant reduction in expression after SCI (log2(fold change) = -1.01) (Figure 3A). High-throughput sequencing revealed significant upregulation (P＜ 0.05), ofcircAbca1after TSCI (Figure 3B).Differentially expressed circRNAs were further screened based on their binding targets. In total, six differentially expressed circRNAs that potentially bindmiR-135b-5pwere screened, includingcircAbca1,circRNA5953,circRNA15386,circRNA3868,circRNA6769, andcircRNA9429(Figure 3C).Based on comprehensive analysis of backsplice junction loci,TargetScan, and miRanda scores, we identifiedcircAbca1as a key circRNA that interacts withmiR-135b-5pto influence pathophysiological processes after TSCI. Reverse-screens identified miRNAs that bind to circAbca1, includingmiR-135b-5p,miR-23a-5p,miR-135a-5p,miR-155-5p, andmiR-138-5p,which also target the bound mRNAs. Together withcircAbca1,these miRNAs/mRNAs were used as the core for constructing the circRNA/miRNA interaction network after TSCI (Figure 3D).

Validation of differentially expressed ncRNAs and confirmation of their relationship

To verify the accuracy of high-throughput sequencing and confirm thecircAbca1/miR-135b-5prelationship, we performed qRT-PCR and fluorescent reporter assays. qRTPCR was used to validate expression ofmiR-135b-5pandcircAbca1, which showed significant differential expression following library sequencing. qRT-PCR results were consistent with our sequencing results (Figure 4AandB). We predicted a targeted binding site forcircAbca1andmiR-135b-5p(Figure 4C), and constructed a luciferase reporter plasmid (Figure 4D) to validate the targeted binding relationship. Luciferase activity of themiR-135b-5pmimic group was significantly lower than that of the negative control group before mutation(P＜ 0.05), but this difference disappeared post-mutation.These results confirm thatcircAbca1can bind tomiR-135b-5p(Figure 4E), consistent with our predicted results.

Validation of functional target gene expression and confirmation of the network relationship

To validate expression of Kruppel-like factor 4 (KLF4) and verify themiR-135b-5p/KLF4relationship, qRT-PCR and fluorescent reporter assays were performed. Subsequently, the potential neuroinhibitory function ofcircAbca1was determined by western blot assay. qRT-PCR was used to validate KLF4 expression, while high-throughput sequencing revealed significant upregulation (P＜ 0.05) ofmiR-135b-5pafter TSCI(Figure 5A), consistent with our qRT-PCR results (Figure 5B).TheLuc-KLF43´UTR-WT plasmid (Figure 5C) was synthesized to validate the targeted binding relationship. The binding site forKLF4andmiR-135b-5pwas predicted by Targetscan software (Figure 5D). Luciferase reporter assay results indicated significantly decreased (P＜ 0.05) luciferase activity in 293T cells co-transfected with MT vector ormiR-135bmimics (Figure 5E).

To confirm the correlation betweencircAbca1,miR-135b-5p, andKLF4in primary mouse dorsal root ganglion spinal neurons, they were transfected with circAbca1 inhibitors.Transfection efficiency and expression ofmiR-135bandKLF4were detected by qRT-PCR. The results revealed that silencing ofcircAbca1increasedmiR-135b-5pexpression, but decreasedKLF4expression (Figure 5F). To further explore the potential neuroprotective functions ofcircAbca1inhibition,a western blot assay was conducted. The results show that GAP43 expression was upregulated aftercircAbca1silencing(Figure 5G).

Figure 1｜Establishment of SCI mouse models and expression profiles of circRNA.

Figure 2｜Enriched GO terms and KEGG pathways of differentially expressed circRNAs in SCI mice.

Figure 3｜circRNA-miRNA-mRNA regulatory interaction network analysis.

Figure 4｜Validation of differential non-coding RNA expression and confirmation of relationships among noncoding RNAs.

Figure 5｜Validation of functional target gene expression and confirmation of the network relationship.

Discussion

TSCI is mainly caused by traffic accidents and falling from high places, but can also be caused by vascular lesions, tumors, and iatrogenic injuries, which often lead to paralysis, unconscious movement, urinary and bowel incontinence, and depression(Cristante et al., 2012). Trauma causes neuronal death and damage to nerve connections, thereby disrupting the upstream and downstream excitation of nerves (Schwab et al., 2018; Song et al., 2018). Physical trauma usually leads to a prolonged secondary injury cascade (Tran et al., 2018) that starts with the expansion of tissue damage from the lesion epicenter; the process is divided into acute, subacute, and chronic phases (Sofroniew, 2018). The subacute phase (from 3 days to 2 weeks post-injury) is considered the best target for future cell-based therapeutic strategies (Rowland et al., 2008).Increasing evidence suggests that non-coding RNAs play an important role in the progression of TSCI (Ning et al., 2014;Kleaveland et al., 2018). Although more than 70% of human genes are transcribed, less than 2% are translated into proteins; most of the remaining genes are transcribed into ncRNAs (Djebali et al., 2012). ncRNAs can be divided into two major categories according to their intrinsic functions:housekeeping ncRNAs, which include nuclear small RNA,ribosomal RNA, and handling RNA; and regulatory ncRNAs,which include miRNAs, circRNAs, and lncDNAs (Qu et al.,2015). Currently, the most comprehensively studied regulatory ncRNA is miRNA, which are widely distributed in varioustissues and body fluids (Ameres and Zamore, 2013; Stoicea et al., 2016). This is particularly obvious in the nervous system,where miRNAs can act as regulatory or even hormone-like factors to affect communication and information transmission between target cells through autocrine or paracrine pathways,which has considerable effects on neurophysiology and neural regeneration (Kleaveland et al., 2018). circRNAs are a new class of ncRNA molecules that play important regulatory roles in cell biology; currently, study of circRNAs is still in an initial stage (Qi et al., 2015; Ren et al., 2019). Unlike other linear structures of ncRNAs, circRNAs exhibit a unique enclosedloop structure and are more stable than linear RNA (Salzman et al., 2013). In the nervous system, circRNAs are expressed in large quantities and highly involved in nerve development,degeneration, and axon growth (Veno et al., 2015; You et al.,2015). Currently, it is believed that the most important role of circRNAs is their activity as miRNA-adsorbing sponges to exert post-transcriptional regulatory functions. For example, recent studies found that thelncCyrano/miR-7interaction network regulates neural activity (Li et al., 2017; Kleaveland et al.,2018; Zhang et al., 2018).

In previous studies, moderate SCI was performed using a modified Allen’s weight drop apparatus (Xie et al., 2018).Excessive damage leads to death or paralysis of mice, but low-intensity damage leads to quick recovery, which reduces observation time. We found that a parameter setting of 6 g/60 mm caused ideal damage. Mice exposed to this condition lose lower extremity motor function within 2 days and gradually recover within 3-21 days after injury (Xie et al., 2018). In the present study, high-throughput sequencing of differentially expressed circRNAs and miRNAs was performed at 5 days post-TSCI. Day 5 post-injury represents the initial period of the subacute phase and scar formation; posttranscriptional regulation may play an important role in pathophysiological processes of the subacute phase. Subsequently, differential expression profiles of circRNAs and miRNAs following subacute injury were obtained.

GO and KEGG analysis of host genes related to differentially expressed circRNAs after TSCI was performed from a macroscopic viewpoint. The results revealed that differentially expressed genes were: 1) mainly involved in transcriptional regulation and protein phosphorylation; 2) abundant in cytoplasm, nucleus, exosomes, and plasma membranes; and 3) play a role in protein, adenosine triphosphate, and metal ion binding. Interestingly, some classical pathways were significantly differentially enriched; for example, pentosephosphate, calcium signaling, p53 signaling, and transforming growth factor beta (TGF-β) signaling pathways. Activation of the pentose-phosphate pathway reduces neuronal oxidative stress (Mashima et al., 2018), while p53 pathway activation drives motor neuron degeneration (Simon et al., 2017) and the calcium signaling pathway regulates synaptic activity (Brini et al., 2014). The TGF-β signaling pathway is one of the most important signaling pathways regulating scar formation and axon regeneration after TSCI, as observed in our previous study (Liu et al., 2018; Wang et al., 2018c).

To determine the specific research core-ncRNA,miR-135b-5pwas identified as the most significantly downregulated miRNA after TSCI.miRNA-135b-5pcan regulate tumor proliferation,migration, and apoptosis in gastric carcinogenesis, cervical cancer, colorectal cancer, mammary carcinomas, and osteosarcoma (Bhinge et al., 2014; Valeri et al., 2014;Chu et al., 2018).miRNA-135b-5palso promotes hypoxiainduced vascular endothelial cell injury (Yang et al., 2018),chondrocyte proliferation (Wang et al., 2018b), pathological cardiac hypertrophy (Chu et al., 2018), and odontoblast-like differentiation (Song et al., 2017). Most importantly, in the nervous system,miR-135b-5pplays a neuroprotective role by targeting GSK3β3 (Duan et al., 2018; van Battum et al.,2018), and can regulate neuroectoderm formation through TGF-β/bone morphogenic protein signaling (Bhinge et al.,2014).miR-135b-5pcan also participate in construction of the network and become its core. It was found that the ncRNA interaction network, made up oflncRNAGAS8-AS1-miR-135b-5p,lncRNAGAS5-miR-135b-5p, andCDR1as-miR-135b-5p, can interfere with cancer progression (Xue et al., 2017; Chen et al., 2019a, b). Our study predicted thatmiR-135b-5pis not just a part of the ncRNA/mRNA interaction network, but at its core. Among all differentially expressed circRNAs capable of binding tomiR-135b-5p,circAbca1exhibited the highest loop index and highest binding rate, and thus was selected for axial prediction.

KLF4is an evolutionarily conserved zinc finger-containing transcription factor that plays key roles in diverse biological processes including proliferation, differentiation, cell growth,and apoptosis (Yuan et al., 2020).KLF4may inhibit axon regeneration and suppress the migration and maturation of palingenetic neurons (Moon et al., 2018). Most importantly,KLF4is the target gene of miR-135b-5p and regulated by the upstream component miR-135b-5p in gastric cancer cells(Chen et al., 2020a, b). Based software prediction, sequencing analysis, and literature retrieval,KLF4was selected as a downstream component in the axial prediction and a complete targetingcircAbca1/miR-135b-5p/KLF4interaction axis was established. Gene expression was verified by PCR,while the results of luciferase reporter assays confirmed that bothKLF4andcircAbca1are target genes of miR-135b.Preliminary verification using circ-inhibitor transfection and detection of functional proteins revealed upregulated expression ofmiR-135b-5pand downregulated expression ofKLF4following circ inhibition, indicating thatcircAbca1could dysregulate the expression of downstream RNAs. In addition to the axial annotation, GAP43 (Li et al., 1996), which plays a critical role in axon growth and regeneration, was also detected. Our results indicated upregulation of GAP43 following transfection with acircAbca1inhibitor, consistent withKLF4downregulation.

In conclusion, the present study preliminarily verified the structure and function of thecircAbca1/miR-135b-5p/KLF4interaction axis, and found thatcircAbca1inhibitors play a neuroprotective role via this interaction network. However,the present study had some limitations. circRNAs commonly bind to several miRNAs; for example, circHIPK3 simultaneously binds to ten miRNAs (miR-654,miR-29b,miR-338,miR-193a,miR-152,miR-193a,miR-379,miR-584,miR-124, andmiR-29a) (Zheng et al., 2016). Moreover, miRNAs bind to several mRNAs (Li et al., 2017; Zhang et al., 2018). Therefore, blocking only circRNA is insufficient to confirm network relationships.Indeed, the regulatory process likely involves a network rather than a single thread, and multiple targeted mRNAs bound by the same miRNA may play synergistic or antagonistic roles in the pathophysiological process after SCI through multiple pathways. Interestingly,miR-135b-5pcan also bind to the TGF-β receptors TGF-βR1 and TGF-βR2 (Li et al., 2015;Wang et al., 2018b). Previously, we confirmed that the TGF-β signaling pathway is important for regulation of scar formation and axon regeneration after TSCI (Liu et al., 2018; Wang et al.,2018c). In subsequent studies, we will focus on the biological roles and mechanisms of thecircAbca1/miR-135b-5p/KLF4interaction network induced by TSCI in this study. Specifically,we will explore and validate the pathological and physiological functions associated with subacute inflammation. We hope that new ncRNA interaction networks and related factors,which can help elucidate the mechanisms underlying TSCI and reveal therapeutic targets and prognostic indicators, can be eventually discovered.

Acknowledgments:The authors express their thanks to Dr. Jia-Nan Chen and Dr. Hong-Fei Li, Shandong University, China and Dr. Yan-Lin Su, Huazhong University of Science and Technology, China for their help in establishing animal models.

Author contributions:Study design: XLH and LL; study implementation and paper preparation: WZW and JL; data analysis: ZDZ, MXL, QL, HXM and HY. All authors approved the final version of this paper for publication.

Conflicts of interest:None declared.

Financial support:This study was supported by the National Natural Science Foundation of China, No. 81874002 (to LL); the Science and Technology Support Project of Sichuan Province of China, Nos. 2018SZ0159 (to LL),2018SZ0246 (to XLH); the Innovation and Entrepreneurship Project of Sichuan Technology Gallery of China, Nos. 2019JDRC0100 (to JL), 2020JDRC0054(to WZW); the National Clinical Research Center for Geriatrics, West China Hospital, Sichuan University, China, Nos. Y2018B22 (to LL), Z20192013 (to JL)and West China Hospital Postdoctoral Research and Development Fund, No.2019HXBH068 (to JL). The funding bodies played no role in the study design, in the collection, analysis and interpretation of data, in the writing of the report,and in the decision to submit the article for publication.

Institutional review board statement:All animal protocols were approved by the Research Ethics Committee of West China Hospital of China (approval No.2017128) on May 16, 2017. The experimental procedure followed the United States National Institutes of Health Guide for the Care and Use of Laboratory Animals (NIH Publication No. 85-23, revised 1996).

Copyright license agreement:The Copyright License Agreement has been signed by all authors before publication.

Data sharing statement:Datasets analyzed during the current study are available from the corresponding author on reasonable request.

Plagiarism check:Checked twice by iThenticate.

Peer review:Externally peer reviewed.

Open access statement:This is an open access journal, and articles are distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike 4.0 License, which allows others to remix, tweak,and build upon the work non-commercially, as long as appropriate credit is given and the new creations are licensed under the identical terms.

Open peer reviewer:Ronak Reshamwala, Griffith University, Australia.

Additional file:Open peer review report 1.