A Potential Competitive Endogenous RNA Pathway Involved in Chronic Spinal Cord Injury.

BACKGROUND Chronic spinal cord injury (CSCI) is a worldwide clinical problem. We aimed to reveal differentially expressed (DE) lncRNAs and to find associated pathways that may function as targets for CSCI therapy. MATERIAL AND METHODS After a CSCI rat model was confirmed by the Basso Beattie Bresnahan (BBB) scale and the Magnetic Resonance Imaging (MRI) test, microarray analysis was used to obtain the expression profile of DE lncRNAs between CSCI rats and corresponding control rats. Then, bioinformatics analyses, including GO and KEGG pathway analysis, DE lncRNAs-mRNAs co-expression analysis, and several databases, were used to examine the function of these DE lncRNAs. Finally, quantitative real-time PCR (qRT-PCR) was used to evaluate the expressions of the 5 most significantly changed lncRNAs, Col6a1, and miR-330-3p. RESULTS Our study identified 1266 DE lncRNAs and 847 DE mRNAs, among which lncRNA6032 was significant up-regulated. Furthermore, the expressions of miR-330-3p and Col6a1 associated with lncRNA6032 were down-regulated and up-regulated, respectively. CONCLUSIONS Our results showed that the abundance of DE lncRNAs may be associated with the risk of CSCI outcome and revealed a potential competitive endogenous RNA (ceRNA) pathway involved in CSCI. Further experiments in vivo and in vitro were essential to uncover the exact mechanism of this ceRNA pathway.

Background

Chronic spinal cord injury (CSCI), a devastating neurological disease, is the chronic phase of spinal cord injury (SCI) resulting in impaired motor, sensory, and autonomic functions. As is well known, SCI includes primary mechanical compression and secondary damage. Our CSCI model is a good simulation of the mechanical compression of the spinal cord and its secondary changes. Because SCI occurs frequently in young adults and dramatically improved treatment therapies have greatly increased life expectancy, the morbidity rate of CSCI is more than 50 times the annual prevalence of acute SCI [1]. In China, CSCI affects millions of people and approximately 120 000 new CSCI cases occur every year [2]. Although a number of studies have documented that inflammation after CSCI is a process contributing to both tissue injury and repair [3,4], the precise cellular and molecular mechanisms, especially the transcriptional or translational regulation of them, is still unknown.

With the rapid development of transcriptomics, we have established that no more than 2% of the human genome encodes for proteins, whereas approximately 75–90% of them are actively transcribed as non-coding RNAs (ncRNAs), most of which are lncRNAs (>200 nucleotides) [5-8]. Recently, increasing research has indicated that lncRNAs have critical roles in all biological processes, including central nervous system development and diseases [9-13]. To date, studies focused on lncRNAs under CSCI are still in their infancy.

Owing to a subdued neural regenerative response, established astrocytic scars, and suspended spontaneous functional recovery, rodents 4 weeks after SCI can be used as a chronic model of SCI [14-16]. In the present study, we conducted a comprehensive transcriptome analysis in rats from the CSCI group and the control (CON) group using Clariom D Rat genechip array and identified differentially expressed (DE) lncRNAs, mRNAs, and pre-miRNAs. In order to discover the function of DE lncRNAs, we then put associated DE mRNAs into the Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene ontology (GO) for further analysis. Then, a co-expression network associated with lncRNAs and mRNAs was constructed. In addition, there is a hypothesis that lncRNAs may function as competing endogenous RNAs (ceRNAs) in modulating the function of mRNAs via competitive binding with miRNAs. The expressions of transcripts of specific mRNAs can be up-regulated through sharing miRNA response elements (MREs) with these lncRNAs. In combination with our previous pre-miRNA profile, we found its corresponding mature miRNA profile through UCSC database and miRBase database. Moreover, we obtained target miRNAs of DE lncRNAs from the UCSC database and RegRNA 2.0 database. Eventually, we discovered that MTA_TC0600001799.mm (which we named lncRNA6032) sequesters miR-330-3p to promote Col6a1 expression, which was established by quantitative real-time PCR (qRT-PCR). In summary, our results reveal a potential ceRNA pathway involved in CSCI.

Material and Methods

Animal

Thirty-two adult male Sprague-Dawley rats (200–250 g) were purchased from Beijing Huafukang Biotechnology Co., Ltd. (Beijing, China). The care of animals and experimental procedures complied with the Animal Ethical and Welfare Committee of the Chinese Academy of Medical Sciences Institute of Radiation Medicine and conformed to internationally recognized guidelines. All rats were maintained under a condition with standard temperature and diurnal light control and had free access to water and food.

Model establishment and motor function

Rats were divided into the CSCI group (n=16) and the CON group (n=16). Simply, rats were intraperitoneally anesthetized using 4% chloral hydrate (10 ml/kg) [17]. Then, rats in the CON group underwent a unilateral laminectomy left at the C6 level. Rats in the CSCI group underwent laminectomy the same as in the CON group and then 1.2×1.2×2 mm dried water-absorbing polyurethane polymer sheets coated with a sustained-release membrane, which was kindly provided by Dr. Yong Hu, and which can reach the maximum volume in 24 h through progressive compression and remain at the maximal expansion, were implanted into the C5 epidural space on the left dorsolateral side. Then, rats were individually housed for 28 days in the environment described previously.

Locomotor activity was evaluated every day after injury using the Basso Beattie Bresnahan (BBB) scale [18]. Double-blinded methods were used by 2 trained researchers and the average score of each group was counted. Moreover, MRI was performed at 28 days after surgery on a 1.5 T clinical MR imaging system (GE Signal HDxt; General Electric Company; USA), with a 4.7-cm surface coil. The ventral recumbent position was preferred for scanning after general anesthesia. The imaging protocol involved T1-weighted imaging (T1WI) (TR/TE range, 480/20 ms; slice thickness 2.0 mm; intersection gap 0.1 mm; FOV 26 cm) and T2-weighted imaging (T2WI) (TR/TE range, 2300/120 ms; slice thickness 2.0 mm; intersection gap 0.1 mm; FOV 26 cm).

Samples collection and total RNAs isolation

Deeply anesthetized rats were perfused through the ascending aorta with phosphate-buffered saline on the 28th day. After that, spinal cords of C5 were harvested. Purified total RNAs were synthesized from spinal cord segments by using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions.

Microarray assay

Synthetization of labelled cDNA and cRNA, Clariom D Rat GeneChips hybridization, and washing, staining, and scanning after hybridization were performed according to Affymetrix protocols. The standardization of hybridization data was processed by Affymetrix GeneChip Operating Software (GCOS). Transcripts with more than 1.1 change folds (FC) and P values of <0.05 were considered as DE. This microarray assay was carried out by Tianjin St. Regis Biotechnology Co., Ltd, Tianjin, China.

Bioinformatics analysis

As the function of lncRNAs was associated with upstream and downstream coding regions (within 10 kb), GO analysis and KEGG pathway analysis were employed for DE mRNAs. Furthermore, co-expression analysis of DE lncRNAs and mRNAs was used if the correlation coefficient between them was more than 0.98. In addition, we obtained DE mature miRNAs according to DE pre-miRNAs through the UCSC database and miRbase database, and target miRNAs of the most significantly changed lncRNAs were obtained from the UCSC database and RegRNA 2.0 database. Among these miRNAs, we searched target mRNAs of common miRNAs by use of several databases, including miRbase and Targetscan.

qRT-PCR validation array

The purified total RNAs described above were reverse-transcribed to cDNAs by using the PrimeScript™ RT Reagent Kit (Takara, Otsu, Japan) in accordance with the standard instructions. Then, we used the DeNovix DS-11 Spectrophotometer (DeNovix, USA) with SYBR Premix Ex Taq™ (Takara, Otsu, Japan) for qRT-PCR testing. The expressions of lncRNAs, including MTA_TC1000000323.mm, lncRNA6032, kakla, lncRNA (chr18: 29647839-29647979) and lncRNA (chr10: 55656404-55656541), miRNA miR-330-3p, and mRNA Col6a1, in the spinal cord were normalized to the expression of the housekeeping gene with the 2−ΔΔCT method [19]. lncRNAs and Col6a were normalized to β-actin and miR-330-3p was normalized to U 6. All experiments were performed 3 times. The primers used are shown in Table 1.

Table 1

List of primers for qRT-PCR.

Gene name	Primer sequence (5′to 3′)	Product length (bp)
MTA_TC1000000323.mm	FW: CCACCTGACGACAATGGAGC	20
	RV: CACAGAAATGGCGTAGGTAGAAG	23
lncRNA6032	FW: TGCGACTTTGTGACTGCTT	19
	RV: TCCTTTCTGGTCCATCTCC	19
kakla	FW: GGGAGGGACTTGAATGTGGA	20
	RV: GCTTGCCTTGGTGAACGCTAT	21
lncRNA (chr18: 29647839-29647979)	FW: GGCCAGCTTTAGCTCAGCG	19
	RV: CGAACCAACCACTCACCAGATAG	23
lncRNA (chr10: 55656404-55656541)	FW: TCCTTACCTGTTCCTCCTCT	20
	RV: TCAGGATGTTGTCAGGTCC	19
miR-330-3p	FW: TGCGGGCAAAGCACAGGGCCTGCAG	25
	RV: CCAGTGCAGGGTCCGAGGT	19
Col6a1	FW: ACTTGATCGTGGTGACTGACG	21
	RV: TGGCGATGATGCTTAGACG	19
β-actin	FW: AGGGAAATCGTGCGTGACAT	20
	RV: CCTCGGGGCATCGGAA	16
U6	FW: GCTCGCTTCGGCAGCACA	18
	RV: AACGCTTCACGAATTTGCGTG	21

FW – forward; RV – reverse.

Statistical analysis

R and Bioconductor programs were utilized to analyze microarray data. The BBB scores were analyzed by multivariate analysis of variance (MANOVA). All continuous variables are expressed as mean ± SD or mean ±SEM, and Mann-Whitney U test and unpaired t test were performed for comparison analysis. The criterion of P<0.05 was statistical significance.

Results

Validation of CSCI model

In our recent research, the animal model was established according to previous reports [20-22]. All rats from the CSCI group were paralyzed in the left forelimb from the first postoperative day, while rats from the CON group were almost normal. As indicated in Figure 1B, locomotor activity of rats in the CSCI group were first worsened and then improved gradually thereafter. However, movements of the CON group rats were regular during the experimental process, as evaluated by BBB scale scores.

Figure 1

(A) T2-weighted transverse (a) cervical MRI of the CON group at C5 spinal cord is normal, while T2-weighted transverse (b) and T1- and T2-weighted sagittal (c, d) cervical MRI of the CSCI group show compression at C5 level with low spinal signal 28 days after surgery. (B) Line graph demonstrates the mean ±SD of the averaged BBB score of the CSCI group (n=16) vs. the CON group (n=16). *** P<0.001.

MRI was used to assess the spinal cord of C5 28 days postoperatively following behavioral tests. Compared to the CON group, the C5 spinal cord of the CSCI group was significantly compressed on the left with weighted sagittal and transverse imaging (Figure 1A).

DE lncRNAs and mRNAs

A graphical outline of the expression characteristics of lncRNAs and mRNAs is shown in heatmap through hierarchical clustering analysis. The heatmap of DE lncRNAs (Figure 2A) or mRNAs (Figure 2B) indicated a high concordance in either CSCI or CON samples. Compared to the control group, 1266 lncRNAs (738 up-regulated and 528 down-regulated) and 847 mRNAs (487 up-regulated and 360 down-regulated), were significantly changed in the CSCI group. We present the top 10 up-regulated and down-regulated lncRNAs and mRNAs in Tables 2 and 3.

Figure 2

Heatmap showing hierarchical clustering of lncRNAs (A) or mRNAs (B) between the CSCI group and the CON group. In clustering analysis, compared to the CON group, up- and down-regulated genes are colored in red and green, respectively.

Table 2

Details of the top 10 up-regulated (A) and 10 down-regulated (B) lncRNAs.

(A) Gene symbol	Identified lncRNA	Gene location	Fold change	P-value
NA	NA	chr10: 55656404-55656541	4.370121079	0.00283592
NA	NA	chr18: 29647839-29647979	2.531533674	3.48684828
MTA_TC1000000323.mm	MTA_TR1000000963.mm	chr20: 44102881-44103700	2.446035343	0.03056308
kakla	kakla.aSep08-unspliced	chr1: 221200474-221204272	2.44155652	0.02868856
MTA_TC0600001799.mm	MTA_TR0600006032.mm	chr4: 29564290-29571817	2.331055367	0.00086609
NA	NA	chr6: 141344474-141345232	2.274962867	0.02572412
XLOC_002175	TCONS_00002938	chr1: 225202604-225214883	2.23071736	0.01884398
MTA_TC0X00000333.mm	MTA_TR0X00000930.mm	chrX: 123705239-123707266	2.082480188	0.00219299
MTA_TC1000001396.mm	MTA_TR1000004367.mm	chr9: 17954218-17954244	2.065710328	0.0063028
MTA_TC0600002866.mm	MTA_TR0600009351.mm	chr4: 133918125-133919091	2.059081008	0.00476291

Table 3

Details on the top 10 up-regulated (A) and 10 down-regulated (B) mRNAs.

(A) Gene symbol	Identified lncRNA	Gene location	Fold change	P-value
Cxcl13	NM_001017496	chr14:15253146-15258220	6.955265872	0.001050302
Postn	NM_001108550	chr2:143656812-143688087	6.076214816	0.000576185
Serpina3n	NM_031531	chr6:128073337-128082355	3.526377745	0.023874686
Ms4a6bl	NM_001006975	chr1:227240417-227254198	2.778223887	0.018778863
Abca1	NM_178095	chr5:69859632-69952466	2.64754434	0.000123318
Ms4a4a	ENSRNOT00000074966	chr1:227217960-227236580	2.557907489	0.043518789
Slfn2	NM_001107031	chr10:70370354-70376718	2.325841943	0.042159484
Hmox1	NM_012580	chr19:14507938-14515457	2.316866342	0.028461287
Ms4a6b	ENSRNOT00000032937	chr1:227798223-227818291	2.227668253	0.030770288
Fbn1	NM_031825	chr3:117569589-117766141	2.076941594	0.04070993

GO and KEGG enrichment analysis

GO and KEGG enrichment analysis of deregulated genes between the CSCI group and the CON group was completed to explore the molecular mechanism of CSCI. GO terms included biological process (BP), cellular component (CC), and molecular function (MF). The top 10 most significantly enriched BP, CC, and MF of up-regulated genes are shown in Figure 3A–3C, respectively. The top 10 most enriched BP, CC, and MF of down-regulated genes are displayed in Figure 3D–3F, respectively. Analogously, the top 10 most significant pathways for up-regulated and down-regulated genes are shown in Figure 4.

Figure 3

Biological functions of DE mRNAs. The significant biological process (A), cellular component (B), and molecular function (C) of up-regulated mRNAs. The significant biological process (D), cellular component (E), and molecular function (F) of down-regulated mRNAs.

Figure 4

KEGG pathway analysis for up-regulated and down-regulated mRNAs. (A) The significant pathway for up-regulated genes. (B) The significant pathway for down-regulated genes.

lncRNA6032-centered lncRNA-mRNA co-expression network

Our previous studies identified several lncRNAs and mRNAs are related to inflammation in CSCI. Therefore, a lncRNA-mRNA co-expression network and a ceRNA network were constructed to investigate potential lncRNAs contributing to inflammation in CSCI rats. By calculating the correlation coefficient value of 519 DE lncRNAs and 562 DE mRNAs, we confirmed that 787 lncRNA-mRNA pairs and 740 lncRNA-mRNA pairs were positively co-expressed and negatively co-expressed, respectively. Based on the 2 distinct co-expression patterns of lncRNA-mRNA associated with inflammation after CSCI, a corresponding lncRNA6032-centered lncRNA-mRNA co-expression network was then constructed (Figure 5).

Figure 5

lncRNA6032-centered lncRNA-mRNA co-expression network. The triangle node and the round node represent lncRNA and mRNA, respectively. The red color shows that the lncRNA or mRNA is significantly up-regulated in the CSCI group compared to the CON group, otherwise indicated as down-regulated with a green color. The positive and negative correlations between lncRNAs and mRNAs are indicated as solid and dashed lines, respectively. The size of the lines represents the significance of the correlations, and the number of single genes that regulate other genes represents the size of the node.

qRT-PCR validation

To validate the reliability of the microarray results, 5 highly expressed DE lncRNAs were selected and validated by qRT-PCR. Col6a1 and miR-330-3p were also chosen to be tested. As shown in Figure 6, qRT-PCR data were in good accordance with the microarray results, further supporting the dependability of the microarray results. Furthermore, compared to the CON group, it was significantly discriminate that the expressions of lncRNA6032 and Col6a1 were considerably higher in the CSCI group while the level of miR-330-3p expression was markedly lower.

Figure 6

Validation of microarray date by qRT-PCR analysis. (A) The 5 highly expressed DE lncRNAs, Col6a1, and miR-330-3p in CSCI in comparison with control by microarray were confirmed by qRT-PCR. The relative expression levels of lncRNA6032 (B), miR-330-3p (C), and Col6a1 (D) in the CSCI group and the CON group by qRT-PCR. Each sample was assessed in 3 technical replicates. The height of each box represents the mean average of sample-specific 2−ΔΔCT values, while associated error bars denote the SEM. * P<0.05, ** P<0.01.

A potential competitive endogenous RNA pathway

Macrophages, which are critical for the inflammatory response, can promote peripheral nerve regeneration and collagen VI generation [23]. Collagen VI, which is essential for nerve myelination, is critical for the migration and polarization of macrophages during peripheral nerve regeneration [24]. In the central nervous system, Collagen VI was also reported to play a neuroprotective role [25,26]. Furthermore, Col6a1, which is expressed by Schwann cells and is an encoding gene for collagen VI, is related to inflammation [27]. According to several databases mentioned above, miR-330-3p is the common target miRNA of lncRNA6032 and Col6a. Generally speaking, we predicted that the suppression of miR-330-3p to Col6a1 can be reduced by lncRNA6032.

Discussion

The pathomechanisms responsible for the high morbidity and disability associated with CSCI are still unknown. lncRNAs, which were initially deemed transcriptional noise [28], were first reported in mice during the extensive sequencing of full-length cDNA libraries [29]. With the development of transcriptomics, it has nevertheless become obvious that lncRNAs regulate gene expression and molecular networks.

Remarkably, 4000–20 000 lncRNA are specifically expressed in the brain [30]. Increasing emerging evidence indicates that lncRNAs exert critical roles in the development of neurological diseases. Therefore, the modulation of lncRNAs may be developed into innovative therapeutic strategies for CSCI. Despite the continuing scientific focus on lncRNAs, very few researchers have concentrated on changes in lncRNAs after CSCI. In our study, a derivation group of 6 rats (the CSCI group and the CON group, n=3 each) was utilized to obtain DE profiles of lncRNAs and mRNAs by using microarray analysis. Then, DE lncRNAs and mRNAs were used for GO and KEGG enrichment analysis and differential co-expression analysis. Furthermore, a ceRNA network was constructed according to bioinformatics. A validation group of 26 rats (the CSCI group and the CON group, n=13) was used to confirm microarray data as well as members of the ceRNA network through qRT-PCR. The design of this experiment supported the robustness of our findings.

The emerging evidence for the function of most lncRNAs was not available to fully understand the regulatory mechanism of CSCI. Consistent with this, the roles of DE lncRNAs identified in the current study were almost unknown. We first used the combination of “GO and KEGG enrichment analysis” with “differential co-expression analysis” [31] as our strategy to research deeper into the biological functions of DE lncRNAs in CSCI. Second, we described 2 different co-expression patterns related to CSCI – positive and negative co-expression networks – which suggested distinct underlying regulation mechanisms. Finally, mRNA Col6a1 associated with inflammation was found to be co-expressed with lncRNA6032. In addition, bioinformatics analysis indicated that miR-330-3p was the target miRNA of both lncRNA6032 and Col6a1. Hence, we speculated that the miR-330-3p-lncRNA6032-Col6a1 interaction may play a critical role in CSCI. Finally, our search of the lncLocator database revealed that lncRNA6032 was located in the cytoplasm.

Increasing reports have proved that lncRNAs can function as ceRNAs [32] or as molecular sponges [33,34] to regulate the function of miRNAs. Our study is the first time to show that lncRNA6032, as a ceRNA, modulates the expression of Col6a1 through competition for miR-330-3p. Based on our findings, significant insights were provided into CSCI pathogenesis, and new treatment strategies could be developed to alleviate it. However, the precise role of these DE lncRNAs in the pathogenesis of CSCI needs to be clarified in our future work. The sample size for the microarray was small, which was a limitation of this study. Research with larger sample sizes and further experiments need to be conducted to confirm this conclusion.

Conclusions

Our results showed that the abundance of DE lncRNAs may be associated with the risk of CSCI outcome and revealed a potential ceRNA pathway involved in CSCI. Further in vivo and in vitro experiments were needed to uncover the exact mechanism of this ceRNA pathway.