Expression profiles of circular RNAs and interaction networks of competing endogenous RNAs in neurogenic bladder of rats following suprasacral spinal cord injury

Experimental animal and SSCI model

Twenty 10–week–old female Wistar rats weighing 200 ± 20 g (Beijing Charles River Laboratories Animal Technology Co., Ltd., Beijing, China) were numbered 1–20 and divided into 4 groups according to a random number table (N =5/group), including the sham group and the SSCI [1–3] groups (1W–, 2W–, and 4W–post–SSCI). All animals were housed in separate cages under well-ventilated conditions with a circadian rhythm of 12/12 h, room temperature of 24 to 26 °C, humidity of 50% to 70%, and free access to water and food. After 1 week of adaptive feeding, all rats were fasted for 24 h and anesthetized with isoflurane inhalation (3%) to reduce the pain or distress of the rats. Rats in the SSCI groups were subjected to Hassan Shaker spinal cord transection (Shaker et al., 2003) to construct the SSCI model. The Institutional Animal Care and Use Committee of Xuanwu Hospital Capital Medical University granted approval for all experimental procedures(No. XW-20210910-1), which were conducted in accordance with the National Institute of Health Guidelines for the Care and Use of Laboratory Animals. The schematic representation of the experimental design was shown in Fig. 1.

The injured segment of this model was uniformly fixed at the T10–T11 segment (Inskip et al., 2009). The method was as follows: The lowest group of floating ribs connected to the 13th thoracic vertebra was used as bony landmarks to locate. Skin preparation and disinfection were performed after determining the T10–T11 level. Longitudinal incision of approximately 2 to three cm was then centered on the vertebral body along the long axis of the spine. The skin and subcutaneous fascia were incised sequentially, and the erector spinae muscles on both sides were bluntly separated to expose the spinous and paraspinal muscles. Muscles were bluntly separated from their attached spinous processes with curved forceps to expose further the spinous processes of T10 and T11 and the adjacent pedicles. The spinal cord was exposed by removing the T10 lamina from the caudal to the cranial side with a microrongeur until the pedicles were on both sides. The spinal cord was hooked out of the intervertebral space with a dental probe and quickly cut off along the hook with a scalpel tip. This can be repeated to ensure complete spinal cord transection. Finally, the two broken ends of the spine were gently lifted with curved forceps to ensure complete spinal cord transection. Injury to the spinal cord capsule was avoided throughout. After observing the heartbeat and respiration of rats without abnormality, the muscles, fascia, and skin were sutured in succession. The incision and surrounding area were disinfected with 5% complexed iodine, and 2 × 10⁵ U of penicillin sodium was injected intraperitoneally. The sham group underwent the same surgical preparation as the experimental group and received laminectomy at the same level. All rats were single-housed after surgery. Clean and dry bedding was changed daily. Reduce the height of feed and water appropriately to ensure that rats can eat independently. Artificial assisted voiding was performed every 8 h in the conscious state of rats using the Crede maneuver until the bladder regained the micturition reflex. Penicillin sodium 2 × 10⁵ U was injected intraperitoneally daily for 5 days after surgery to prevent infection.

Inclusion and exclusion criteria

Inclusion criteria: 1. Hind limb motor function: rats’ hind limbs in the experimental group were observed if they participated in walking. If the hindlimb was only dragged while walking on the forelimb and could not participate in walking, the Basso, Beattie, and Bresnahan (BBB) score for SSCI in rats was 0. It was then considered a successful SSCI model. 2. Bladder voiding function: In the surgical anesthesia and spinal shock period, the bladder was observed to see if the state of retention and bladder distension were obvious, and spontaneous urination could not be performed. However, the bladder capacity decreased after the spinal shock period, and its distension was not obvious. No manipulation was required to assist urination, and the lower abdomen and bedding material in the cage were humid, indicating that the detrusor muscle of the bladder in rats had produced uninhibited contraction.

Exclusion criteria: Model rats were not included in the experiment if they showed the voluntary movement of both hind limbs, spontaneous urination, autophagy, or death.

Bladder tissue collection

For ethical reasons, all animals were euthanized according to AVMA guidelines. We administered intraperitoneal injections of pentobarbital sodium at a dosage of 150 mg/kg to perform euthanasia, thereby reducing the pain or distress of the rats. Bladder tissues were collected at 1, 2, and 4 weeks after modeling in the SSCI [1–3] groups, respectively, and at 4 weeks after surgery in the sham group. Each tissue was divided into three equal fractions. One-third of the samples were preserved in 4% paraformaldehyde, embedded in paraffin, and then sectioned for morphological observation using hematoxylin-eosin (HE) and Masson staining. Another one-third was used to extract total RNA for transcriptome high-throughput sequencing. The last one-third was used for quantitative real-time polymerase chain reaction (qRT-PCR) validation.

HE staining and Masson staining

Pathological analyses were performed to assess the size of smooth muscle cells in the bladder and the presence of bladder fibrosis. The bladder tissues were fixed overnight with a 4% paraformaldehyde solution, then dehydrated in a series of ethanol and cleared in xylene. The tissues were subsequently embedded in paraffin and cut into sections with a thickness of 6 micrometers. HE staining detected the nuclei, bladder smooth muscle cells, and bladder morphology. Masson staining showed collagen fibers (blue) and muscle tissue (red) to detect the severity of bladder fibrosis. The staining signals were quantitatively analyzed using Image J software (version 1.53q).

RNA isolation and cDNA library construction and Illumina sequencing

Three bladders were obtained from each group for whole transcriptome sequencing. TRIzol reagent (Invitrogen, Carlsbad, CA, USA) was used to extract total RNA from each tissue. RNA concentration in each tissue was measured using NanoDrop ND-2000 (NanoDrop, Wilmington, DE, USA). The RNA Integrity Number (RIN) was determined using the Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA, USA); all tissues had RIN ≥7.0. Twelve cDNA libraries were generated, each containing three different tissue samples to analyze the expression patterns of both mRNAs and circRNAs. According to the protocol, ribosomal RNAs (rRNAs) were depleted using TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, CA, USA). Afterward, the TruSeq Stranded Total RNA Library Prep Kit (Illumina, San Diego, CA, USA) was utilized to construct the cDNA libraries. Quantification of libraries was determined using the Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA, USA). Each group had 3 biological replicates when performing RNA sequencing (RNA-seq) at approximately 72X sequencing depth. The statistical power of this experiment was 0.86 as calculated in RNASeqPower (https://doi.org/doi:10.18129/B9.bioc.RNASeqPower). The Illumina Hiseq 4000 Sequencer (Illumina, Shanghai, China) was used to conduct paired-end reads mode sequencing of 150 bp on the final cDNA libraries.

Discovery of the mRNAs and circRNAs

The raw reads underwent quality control using Trimmomatic software (version 0.36) with the default quality threshold to remove the containing adaptors, excessive N content of unknown bases, and low-quality tags. After this step, the clean reads were ready for mapping and subsequent analysis. FastQC (version 0.11.5) was used to establish data quality sequence that persisted after filtration. High-quality clean reads were mapped to the reference genome Rnor_6.0 by Hisat2 software (version 2.0.5) with the rna-strandness parameter. The mapped reads from each library were assembled using StringTie2 software (version 1.3.3b), and a unified transcriptome was created by combining all the tissues. CIRI2 software (version 2.0.3) and circAtlas (version 2.0) were used to identify the circRNAs.

Functional analysis of DEMs and DECs

The readings of the SSCI [1–3] groups were individually compared to the sham group. DESeq2 software (version 1.18.0) was used, which employs the negative binomial distribution to determine the differential expression between the SSCI and sham tissues. Volcano plots and hierarchical clustering were used to characterize and identify the DEMs and DECs by R software (version 3.22.3) with absolute fold change (FC) >2 (p < 0.05). The logarithmic transformation was applied to the FC values of the matched groups. A significance level of p < 0.05 was used as the threshold for identifying DEMs and DECs. The RNA expression of the sham group was used as the denominator during FC computation. The DEMs and DECs with log2 (FC) <−1 were considered down-regulated, while those with log2 (FC) >1 were considered up-regulated. The DEMs and DECs associated with each NB stage were subjected to Venn diagram analysis. Thus, we determined the commonly expressed DEMs and DECs that potentially up-regulated or down-regulated gene expression in the three-time points of NB.

The Miranda software (version 3.3a) predicted the DECs-targeted miRNAs to determine the functions of DECs in bladder fibrosis and NB following SSCI. The target genes of DECs were analyzed using R software (package clusterProfiler) to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. The target miRNAs of DECs and the target mRNAs of miRNAs were predicted using TargetScan, MiRanda, and RNAhybrid. A Cytoscape-based approach (http://cytoscape.org/) was used to construct the circRNA–miRNA–mRNA interaction network.

qRT-PCR Validation

The top 10 co-expressed DECs in each NB stage were subjected to qRT-PCR to verify their expression. TRIzol reagent was used to extract total RNA from bladder tissues. After RNA extraction, rRNAs were depleted using TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA). Subsequently, the RNA samples were fragmented and reverse transcribed into cDNA using random primers via the PrimeScript RT reagent kit(Takara Bio Inc., Japan). Then, qRT-PCR was performed using TB Green Premix Ex Taq II (Takara Bio Inc., Shiga, Japan). Gapdh served as the internal control for normalization. All primers for qRT-PCR were synthesized by Tianyi Huiyuan Co. (Beijing, China).

Statistical analysis

Statistical analyses were performed using SPSS software (version 24.0) and GraphPad Prism (version 8.0). The mean ± standard error of mean(SEM) was used to present normally distributed data. The t-test was utilized to determine significant differences in circRNAs and mRNAs expressions between groups, and statistical significance was set at a p-value of <0.05.

Validation of NB following SSCI model in rats

HE and Masson staining (Figs. 2A–2H) were used to verify whether the SSCI rat model was successfully established and evaluate the bladder fibrosis degree. Furthermore, we conducted a detailed quantitative analysis of the staining signals using Image J software (Fig. 2I). The results showed that only a small amount of collagen fibers was seen in the muscular layers of the bladder in the sham group. In the SSCI [1–3] groups, the normal fibrous connective tissue in the bladder wall gradually disordered; the layered structure of the bladder wall gradually disintegrated, and the bladder wall first thinned and then thickened with the time of injury; the lamina propria decreased, smooth muscle hypertrophy and thickening; intermuscular collagen fibers gradually increased, and the staining gradually deepened. According to all the above features, NB following the SSCI model had been established successfully.

Results of sequencing and characteristics of transcripts

Using the Illumina platform, 12 cDNA libraries have been sequenced, resulting in approximately 187 Gb of sequenced data. The clean reads of each tissue were compared with the Rnor_6.0 reference genome, and the alignment rate ranged from 93.54 to 94.76% (Table S1). The valid ratio for each library is displayed in Table S2. An examination was conducted on mRNAs and circRNAs expression profiles in both the sham and SSCI [1–3] groups. The results showed that 22,601 mRNAs (Fig. S3) and 7,250 circRNAs (Fig. S4) were identified in all bladder tissues.

The mRNAs expression profiles in sham group and SSCI [1–3] groups

Volcano maps were used to assess the locations of differential mRNAs of the SSCI [1–3] groups vs. the sham group (Figs. 3A–3C). Compared with the sham group, 3,255 (1,885 up-regulated and 1,370 down-regulated), 3449(2052 up-regulated and 1,397 down-regulated), and 884 (254 up-regulated and 630 down-regulated) DEMs were detected in the SSCI–1, SSCI–2, and SSCI–3 groups, respectively. The SSCI [2–3] groups also exhibited some DEMs compared to the SSCI–1 group. Specifically, there were 1,056 and 3,228 DEMs in the SSCI–2 and SSCI–3 groups, with 636 and 1,077 DEMs showing up-regulation, and 420 and 2151 DEMs showing down-regulation, respectively. Compared with the SSCI–2 group, the SSCI-3 group had 3,024 DEMs, with 937 showing up-regulation and 2,087 showing down-regulation (Fig. 3D). Clustered heatmap showed differences in mRNAs expression patterns between the SSCI [1–3] groups and the sham group (Fig. 3E). Venn diagrams of the up-regulated and down-regulated DEMs showed that 420 (115 up-regulated, 305 down-regulated) coexisted in the intersections of the SSCI [1–3] groups vs. the sham group (Figs. 3F–3G). The specific comparisons of DEMs between each group could be found in Fig. S3.

The circRNAs expression profiles in sham group and SSCI [1–3] groups

Volcano maps were used to assess differential circRNAs locations of the SSCI [1–3] groups vs. the sham group (Figs. 4A–4C). The circRNAs expression profiling data showed that 1,339, 1,324, and 1,151 DECs were differentially expressed between the SSCI [1–3] groups and the sham group, of which 408, 349, and 532 were up-regulated, and 931, 975, and 619 were down-regulated, respectively. The SSCI [2–3] groups also exhibited some DECs compared to the SSCI–1 group. Specifically, there were 746 and 1,135 DEMs in the SSCI–2 and SSCI–3 groups, with 332 and 751 DECs showing up-regulation, and 414 and 384 DECs showing down-regulation, respectively. Compared with the SSCI–2 group, the SSCI–3 group had 1,120 DECs, of which 816 were up-regulated, and 304 were down-regulated (Fig. 4D). According to the clustered heatmap, circRNAs expression patterns varied between the SSCI [1–3] groups and the sham group (Fig. 4E). Venn diagrams of the up-regulated and down-regulated DECs showed that 412 (52 up-regulated, 360 down-regulated) coexisted in the SSCI [1–3] groups vs. the sham group intersections (Figs. 4F–4G). Additional analyses of the DECs were conducted to investigate the molecular features of circRNAs, focusing on their length and chromosomal distribution. Fig. 4H revealed that DECs majority were derived from chromosome 1 (10.83%; 785/7250), followed by chromosome 2 (8.0%; 580/7250), and chromosome 5 (6.94%; 503/7250), in descending order. Furthermore, Fig. 4I demonstrated that a significant proportion of the DECs had lengths of less than 2,000 nucleotides (85.94%; 6233/7250). Tables 1 and 2 presented the top 10 up-regulated and down-regulated DECs, respectively, that are co-expressed in all SSCI [1–3] groups when compared with the sham group. The specific comparisons of DECs between each group could be found in Fig. S4.

Validation of the expression of circRNAs by qRT-PCR

The RNA-seq analysis showed that circRNA3621, circRNA5775, circRNA0617, circRNA2378, and circRNA0613 were the five most significantly up-regulated DECs common to SSCI [1–3] groups compared with the sham group. Moreover, circRNA3111, circRNA0586, circRNA1617, circRNA4426, and circRNA5665 were the five most significantly down-regulated DECs common to SSCI [1–3] groups compared with the sham group. The above 10 circRNAs were selected for qRT-PCR to confirm the RNA-seq results’ accuracy, and the sequences were listed in Table 3. The qRT-PCR experiment results, presented in Fig. 5, were consistent with the RNA-seq data, thus affirming the RNA-seq results precisions.

The feature and the potential function of DEMs and DECs

GO enrichment analysis was performed, during which biological process (BP), cellular component (CC), and molecular function (MF) were described to filter the key regulators and pathways of the DEMs and the hosting genes of the DECs in the process of NB following SSCI.

For DEMs, the SSCI [1–3] groups were compared to the sham group (Fig. S6); the top 10 up-regulated and the top 10 down-regulated GO terms of DEMs were shown in Figs. 6A–6F. In the up-regulated GO terms, the BP category was mostly enriched in “inflammatory response”, “inflammatory response”, and “chronic inflammatory response”. The CC category was mostly enriched in “extracellular space”, “extracellular space”, and “extracellular space”; the MF category was mainly involved in “integrin binding”, “cytokine activity”, and “serine-type endopeptidase activity”. In the down-regulated GO terms, the BP category was mostly enriched in “potassium ion transmembrane transport”, “potassium ion transport”, and “positive regulation of synapse assembly”; the CC category was mostly enriched in “voltage-gated potassium channel complex”, “integral component of the plasma membrane”, and “anchored component of membrane”; the MF category was mostly involved in “calcium ion binding”, “calcium ion binding”, and “calcium ion binding”. In addition, comparing the SSCI [1–3] groups to the sham group (Fig. S7), the top 20 up-regulated and down-regulated KEGG pathways are shown in Figs. 7A–7F. Among those, the most significantly up-regulated pathways were “Cytokine-cytokine receptor interaction”, “Cytokine-cytokine receptor interaction”, and “Cytokine-cytokine receptor interaction”, respectively. The most obvious down-regulated pathways were the “Neuroactive ligand–receptor interaction”, “Calcium signaling pathway”, and “Neuroactive ligand–receptor interaction”, respectively.

The top 30 up-regulated and the top 30 down-regulated GO terms of DECs are shown in Figs. 8A–8F. According to GO analysis (Fig. S8), the terms that contained the most DECs were “negative regulation of neuron migration”, “negative regulation of cell migration”, and “negative regulation of receptor biosynthetic process”, respectively, in the up-regulated BP category. The terms that contained the most DECs were “extracellular region”, “trans-Golgi network membrane”, and “postsynaptic density”, respectively, for up-regulated CC category. For up-regulated MF, the terms that contained the most DECs were “iron ion binding”, “pre-mRNA intronic binding”, and “RNA polymerase II transcription factor activity, sequence-specific transcription regulatory region DNA binding”, respectively. In the down-regulated GO terms, the BP category was mostly enriched in “adenylate cyclase-activating G-protein coupled receptor signaling pathway”, “cAMP biosynthetic process”, and “intracellular signal transduction”; the CC category was mostly enriched in “endocytic vesicle”, “guanylate cyclase complex, soluble”, and “sarcolemma”; the MF category was mostly involved in “calcium- and calmodulin-responsive adenylate cyclase activity”, “calcium- and calmodulin-responsive adenylate cyclase activity”, and “calcium- and calmodulin-responsive adenylate cyclase activity”. In addition, comparing the SSCI [1–3] groups to the sham group (Fig. S9), the top 20 up-regulated and the top 20 down-regulated KEGG pathways are shown in Figs. 9A–9F. Among those, the most significantly up-regulated pathways were involved in “Vasopressin-regulated water reabsorption”, “Notch signaling pathway”, and “Glutamatergic synapse”, respectively. The most obvious down-regulated pathways were involved in the “Phospholipase D signaling pathway”, “Phospholipase D signaling pathway”, and “Bile secretion”, respectively. Go and KEGG analysis indicated that an increasing quantity of cells with immune function compared to the sham group, and the Notch signaling pathway might be the main cause of NB following SSCI.

The circRNA–miRNA–mRNA network construction

Based on comparing the SSCI groups and the sham group at three different time points, we constructed the circRNA–miRNA–mRNA co-expression network to reveal the complex interactions among circRNAs, miRNAs, and mRNAs. Those DECs predicted miRNAs interacting with them, and the ceRNAs network was calculated and visualized according to those miRNAs (Fig. S10). These circRNAs and their neighbors were found to form a complex module when they were mapped onto a co-expression network. As shown in Fig. 10, circRNA1249 and circRNA2239 were the network’s largest node, indicating that the circRNA1249-mediated and circRNA2239-mediated network may play a crucial role in the NB process following SSCI.