LncRNA/miRNA/mRNA ceRNA network analysis in spinal cord injury rat with physical exercise therapy

Establishment of spinal cord injury models

We obtained all the adult female Sprague Dawley rats (220–250 g) from the Animal Care Facility of Nantong Medical University (Nantong, China). All animal experiments were performed in accordance with institutional norms and were permitted by Nanjing Medical University’s Animal Care and Use Committee (S20200317-024). All animals were bred and housed (five per cage) in a controlled environment with an ambient temperature of 23 ± 2 °C and free access to water and food was given to animals for a 12-hour light-dark cycle. Rats received a hemisection SCI at the spinal level of Thoracic 10 as described in our previous report. The hemisection surgery was carried out by an expert in the SCI rat model to eliminate the error caused by the incision. Following surgery, an intramuscular injection of the Penicillin 20IU/d (i.m.) was given to the animals for 7 days after the operation. Rats were given bladder care three times each day for two weeks, or until the recovery of bladder control. Baytril (0.06 mg/kg) was injected subcutaneously two times per day for 7 days if there was any sign of infection. Certain criteria were considered for animals being included in or excluded from the experiment. The animal that met any two of the group I criteria ((a) rough coat and unkempt, (b) eyes completely or partially closed for 10 minutes, (c) markedly diminished resistance to being handled (grimace response), (d) markedly decreased movement/lethargy, (e) hunched posture, and (f) distended abdomen), were excluded and euthanized. Euthanasia was performed on mice who met one of the group II criteria ((a) inability to eat or drink, (b) moribund/unresponsive, (c) failure to right itself when placed on its back, (d) dyspnea, or (e) 15% or greater reduction in body weight). The rats were sacrificed two weeks after the second operation and the epicenter of the T10 spinal cord was collected for sequencing and q-PCR analysis. All mice were euthanized humanely, via carbon dioxide asphyxiation (all rats were included in the analysis).

Exercise therapy administration

Animals were given body weight-supported Treadmill Training (TMT) for two weeks starting one week after SCI. The feet of the rats were secured to the pedals and suspended horizontally. The treadmill speed was kept at 6 m/min, and each exercise bout comprised two 20-minute workout intervals per day, five days per week for two weeks. The weight-supported range throughout the training was set at 20%–40% of the rat’s weight, depending on the rat’s functional state. Treatments were administered daily in the same arrangement and at the same time, and cages locations were kept at the same positions throughout the experiment to minimize confounders.

Sample size for RNAseq studies

The RNA-seq study’s error is caused by two factors: the sequencing’s technical variability , and the biological variability of the specimens being compared within groups. For human samples, this value is frequently 4-1, but for inbred animal lines, it is frequently 1 or less. That is, a 2-fold difference is observed in mean expression between treated and untreated samples, while also seeing a 50% variation in expression within the control or treatment groups. This would correspond to an effect size of 2 and a biological coefficient of variation in group 1 /or both (CV) of 0.5. In general, sequencing depths greater than 5/CV² result in only minor increases in study efficiency and/or power, but adding more samples is always efficient.

Depth is a required argument; any one of the others may be left missing and the function will solve for it. An argument may be a vector, in case a vector of values is returned. If multiple arguments are vectors a matrix or array of results is returned. Values for alpha and power would be 0.05 and 0.9. The effect parameter specifies biological effect that, if it were true, the experimenter would want to be able to detect; commonly used. The statements that A has twice the expression of B has half of A are symmetric, values of 1.5 and 2 will the same result. The estimated CV of expression within group may be the most difficult parameter to choose the in depth discussion recommended values for types of data. Samples sizes in the two groups are equal (n2 = n) if a second sample is not . Likewise the coefficients of variation in the two groups are assumed to be equal if cv2 is not specified. We show the tables of sample size results at different depths according to each power value through https://git.bioconductor.org/packages/RNASeqPower (Table S1–S3).

miRNA extraction and sequencing

In total, six samples (three from the SCI group and three from the SCI plus TMT group) were selected for high-throughput RNA sequencing. For each sample, three random rats spinal cord tissues in the same group were mixed together before RNA extraction. A five mm spinal cord was taken from the intact and injured site, the spinal cord segments were rapidly frozen on dry ice following the storage at −80 °C before RNA extraction. After total RNA was extracted by Trizol reagent kit (Invitrogen, Carlsbad, CA, USA), the RNA molecules in a size range of 18–30 nt were enriched by polyacrylamide gel electrophoresis(PAGE). After adding the 3′ adapters, the 36–44 nt RNAs were enriched. Additionally, the 5′ adapters were ligated to the RNAs. The ligation products were reverse transcribed using PCR amplification. The PCR products, which ranged in size from 140 to 160 bp, were enriched to create a cDNA library and sequenced on an Illumina HiSeqTM 2500 (Gene Denovo Biotechnology Co, Guangzhou, China). Results obtained from sequencing contained noise which would affect the following assembly and analysis. Thus, to get pure tags, raw reads were further filtered according to the fastqc software. Then, all purified tags were aligned with GeneBank, Rfan, and Genome. miRNA including existing miRNA, known miRNA, and novel miRNA was identified by software Mireap_v0.2. After tags were annotated as mentioned previously, the annotation results were determined in this priority order: rRNA etc > exist miRNA > exist miRNA edit > known miRNA > repeat > exon > novel miRNA > intron. The tags that cannot be annotated as any of the above molecules were recorded as unann. Correlation analysis of two parallel experiments provides the evaluation of the reliability of experimental results as well as operational stability. The correlation coefficient between two replicas was calculated to evaluate repeatability between samples. The closer the correlation coefficient gets to 1, the better the repeatability between two parallel experiments. We identified miRNAs with a fold change ≥2 and P value <0.05 in a comparison as significant DE miRNAs.

mRNA and LncRNA bioinformatics analysis

Total RNA extraction was performed according to the manufacturer’s procedure using the Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) . The integrity of the RNA was determined using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and validated using RNase free agarose gel electrophoresis. After extracting total RNA, rRNAs were removed to retain mRNAs and ncRNAs. The enriched mRNAs and ncRNAs were fragmented into short fragments by using fragmentation buffer and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP (dUTP instead of dTTP) and buffer. Next, the cDNA fragments were purified with QiaQuick PCR extraction kit (Qiagen, Venlo, the Netherlands), end repaired, poly(A) added, and ligated to Illumina sequencing adapters. Then UNG (Uracil-N-Glycosylase) was used to digest the second-strand cDNA. The digested products were selected on the bases of size by agarose gel electrophoresis, PCR product, and sequenced using Illumina HiSeqTM 4000 (or other platforms) by Gene Denovo Biotechnology Co. (Guangzhou, China). Clean reads were acquired using fastqc software (version 0.18.0) and reads were then mapped to the ribosomal RNA (rRNA) database using the short read alignment tool Bowtie2 (version 2.2.8). After that, the reads that had been mapped to rRNA were deleted. The remaining reads were further used in the assembly and analysis of the transcriptome. An index of the reference genome was built, and paired-end clean reads were mapped to the reference genome using HISAT2 (version 2.1.0). The reconstruction of transcripts was carried out with software Stringtie. To identify the new transcripts, all reconstructed transcripts were aligned to reference genome and were divided into twelve categories using Cuffcompare. Transcripts with one of the class codes “u,i,j,x,c,e,o” were defined as novel transcripts. Then we used the length of the transcript which was longer than 200bp and the exon number was more than 2 to identify reliable novel genes. Two software i.e., CNCI (version 2) and CPC (version 0.9-r2) (http://cpc.cbi.pku.edu.cn/), with default parameters, were used to assess the protein-coding potential of novel transcripts. The intersection of both non protein-coding potential results was chosen as long non-coding RNAs.

Transcripts abundances were quantified by software StringTie is a reference-based approach. For each transcription region, an FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify its expression abundance and variations. The coding and non-coding RNA transcripts that were differently expressed were investigated, accordingly. DESeq2 was then used to compare the differential expression between two groups (or by edgeR between two samples). The genes/transcripts with the parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≥2 were considered as differentially expressed genes/transcripts.

Real-time qRT-PCR validation

For expressions quantification, total RNA from spinal cord samples was extracted by using TRIzol reagent (Pufei biological, USA). Extracted RNA was converted into cDNA using a reverse transcription kit (Vazyme, Nanjing, China) as per the instructions of the manufacturer. Internal control for the relative expressions of the gene was set as glyceraldehyde-3-phosphate dehydrogenase gene (GAPDH). Relative levels of gene expression were calculated using the 2-ΔΔCt method (n = 3). The primers used are given in Table 1.

Constructing the CeRNA network

For evaluation DELs-DEMs interactions the starBase database (version 2.0; https://starbase.sysu.edu.cn/starbase2/) was used (Li et al., 2014), integrating with miRNA-mRNA interactions to create networking of DEL-DEM-DEG ceRNA employing Cytoscape software (version 3.6.1; http://www.cytoscape.org). Similarly, the DECs and DEMs interactions were coupled and then combined with the miRNA-mRNA interactions using Cytoscape software for creating the network of DEC-DEM-DEG ceRNA.

Functional enrichment and annotation

Gene ontology (GO) analysis is a popular bioinformatics tool for annotating and analyzing transcriptome data. KEGG is a database that helps researchers better understand the biological system’s advanced functions and applications. The data acquired from GO and KEGG pathway database analysis for DEGs were processed by using the Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.ncifcrf.gov/). For the genome study, REACTOME (a pathway database) was employed to illustrate the molecular pathway information.

Statistical analysis

To conduct statistical analysis, we used the GraphPad Prism software. The means ± standard deviations (mean ± SEM) were used to express all of the data. For RT-PCR, three technical replicates were collected from each independent experiment to account for variability in the sample. The significances among multiple groups were determined by One-way Analysis of Variance (ANOVA) followed by Bonferroni post hoc test. P value < 0.05 was considered statistically significant.

The suppression of the post-SCI inflammatory response by physical exercise therapy

Rats were divided into three groups: (i) Sham; (ii) SCI; (iii) SCI plus TMT (TMT). Firstly, we successfully constructed an SCI model (Fig. S1). In previous studies, SCI-induced microglial activation and successive release of inflammatory factors i.e., tumor necrosis factor (TNF), interferon (INF), and interleukin (IL) have been shown to cause direct death of neurons however induction in vascular endothelial cells for the expression of a range of chemotaxis and cell adhesion molecules (Li et al., 2020a; Li et al., 2020b; Li et al., 2020c; Orr & Gensel, 2018). In our study, within three weeks of post-SCI, the levels of mRNA expressions of IL-6, TNF-α, and IL-1β were increased. Treadmill training was performed after one week of SCI to assess the anti-inflammation potential of physical exercise therapy. Two weeks later, the levels of mRNA of the three pro-inflammatory cytokines at the site of the lesion were examined (Fig. 1A). All three cytokines were markedly suppressed after treatment with physical activity (Figs. 1B–1E).

Differential expression analysis

As per the pre-set threshold (FDR value <0.05 and |FC| ≥ 2), we have identified a whole of 76 differential lncRNAs between SCI as well as in control samples, comprising 41 upregulated and 35 downregulated lncRNAs. Whereas 33 differential miRNAs were identified, containing 21 up-regulated and 12 downregulated miRNAs. In the account of mRNAs, 19 differential mRNAs were generated, comprising 11 upregulated and eight downregulated genes. All differential LncRNAs, miRNAs, and mRNAs are presented in Fig. 2 and Tables 2–4.

GO ontology

Gene ontology was used to determine which molecular functions and biological processes were overexpressed in DEGs in different analyses employing DAVID 6.8. The top 20 GO terms were shown in Table 4. biological process and 12 cellular components were significantly (p < 0.05) identified in “spinal cord injury rats vs Control”. Most of these processes were related to cellular components, for example, Intracellular organelle, whole membrane, synapse, neuron projection, mitochondrion, macromolecular complex, organelle membrane (Table 5). The DEGs are mainly involved in eight molecular functions, 13 cellular components, 21 biological processes (Figs. 3–4). In the pathophysiology of SCI, these functions are critical for the control of differently expressed mRNAs.

KEGG pathway enrichment analysis of the DEGs

The pathways related to DEGs in SCI rats were mostly enriched, according to KEGG analysis. The differentially identified pathways include metabolic pathways, proteoglycans in cancer, phagosomes, cell adhesion molecules, fluid shear stress, and atherosclerosis, ubiquitin mediated proteolysis, cardiac muscle contraction, ECM-receptor interaction, p53 signaling pathway, lysine degradation, SNARE interactions in vesicular transport, circadian rhythm, and ribosome. The fold enrichment of DEGs involved in each pathway was shown by a bar chart(Fig. 5 and Tables 6–7). According to reactome pathway enrichment analyses, PSGs bind proteoglycans, and TGF-beta1 was the most significantly affected phase in SCI. This was supported by the findings of the GO enrichment analysis.

ceRNA network

Nine DEMs were projected to control 45 DELs using the starBase database, and this information was applied for the development of the lncRNA-miRNA-mRNA ceRNA network, which was then interconnected to the miRNA-mRNA network. In this network, there were 73 nodes (20 DEMs; 45 DELs; 8 DEGs), and 246 interactions (21 DEL-DEM and 225 DEM-DEG interactions) (Fig. 6). Particularly, endosome-related signaling pathways were considerably enriched in the genes functional analysis in the network of lncRNA-associated ceRNA (Table 8). The detail of ceRNA network show in Table S4.