Chronic constriction injury of sciatic nerve changes circular RNA expression in rat spinal dorsal horn

Abstract

Background

Mechanisms of neuropathic pain are still largely unknown. Molecular changes in spinal dorsal horn may contribute to the initiation and development of neuropathic pain. Circular RNAs (circRNAs) have been identified as microRNA sponges and involved in various biological processes, but whether their expression profile changes in neuropathic pain condition is not reported.

Methods

To test whether neuropathic pain influences circRNA expression, we developed a sciatic chronic constriction injury (CCI) model in rats. The CCI ipsilateral spinal dorsal horns of lumbar enlargement segments (L3–L5) were collected, and the total RNA was extracted and subjected to Arraystar Rat circRNA Microarray. Quantitative real-time polymerase chain reaction (qPCR) was used to confirm the circRNA expression profile. To estimate functions of differential circRNAs, bioinformatics analyses including gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes Pathway analyses were performed for the top 100 circRNAs and circRNA–microRNA networks were constructed for the top 10 circRNAs.

Results

circRNA microarrays showed that 469 circRNAs were differentially expressed between CCI and sham-operated rats (fold change ≥2). In all, 363 of them were significantly upregulated, and the other 106 were downregulated in the CCI group. Three of them (circRNA_013779, circRNA_008008, and circRNA_003724) overexpressed >10 times after CCI insult. Expression levels of eight circRNAs were verified using qPCR. GO analysis revealed that thousands of predicted target genes were involved in the biological processes, cellular component, and molecular function; in addition, dozens of these genes were enriched in the Hippo signaling pathway, MAPK signaling pathway, and so on. Competing endogenous RNAs analysis showed that circRNA_008008 and circRNA_013779 are the two largest nodes in the circRNA–microRNA interaction network of the top 10 circRNAs.

Conclusion

CCI resulted in a comprehensive expression profile of circRNAs in the spinal dorsal horn in rats. CircRNAs in the dorsal horn could be helpful to reveal molecular mechanisms of neuropathic pain.

Keywords:

• neuropathic pain
• circular RNA
• circRNA–microRNA interaction
• microarray

Introduction

Millions of people are suffering from neuropathic pain worldwide, which is still a major health and economic burden.^1,2 The exact molecular mechanisms of neuropathic pain are largely unknown and the elucidation of its molecular mechanisms will be beneficial for the development of mechanism-oriented treatment.³

The spinal dorsal horn harbors secondary neurons that receive nociceptive stimuli. It is reported that the dorsal horn plays an important role in the generation of neuropathic pain and the occurrence of central sensitization of pain.^4–7 Dorsal horn is also the location where inhibition of the transmission of nociceptive stimuli takes place because the descending inhibitory neuronal terminals end here;^8,9 therefore, it is presumed that regulatory molecules play a role in the dorsal horn of the spinal cord in neuropathic pain, because neuroplastic changes in the peripheral (such as the dorsal root ganglion) and central nervous system (such as the spinal dorsal horn and brain) of the pain processing network take part in the development of pathological pain.^10,11

MicroRNAs are noncoding posttranscriptional molecules and important biomarkers, which modulate numerous pathophysiological processes including neuropathic pain.¹² For example, microRNA-7a,¹³ microRNA-1,^14,15 and microRNA-30b¹⁶ in the dorsal root ganglion have been shown to be involved in the induction of neuropathic pain. Although one study reported that chronic constriction injury (CCI) did not lead to relevant changes in spinal miRNA expression levels in rats, Genda et al¹⁷ detected >100 differentially expressed microRNAs in the dorsal horn of CCI rats. Another study reported that microRNA-203 in the dorsal horn decreased as many as 10 fold in CCI rats.¹⁸

Circular RNAs (circRNAs) are another class of noncoding RNAs which interact with microRNAs.¹⁹ CircRNAs act as miRNA sponges^20–22 and regulate gene expression through a circRNA–microRNA–mRNA pathway.^23,24 Accordingly, we believed that altered expression of circRNAs may accompany the development and progress of neuropathic pain. To date, studies have been performed probing the promising role of circRNAs as biomarkers of diseases such as Alzheimer’s disease,²⁵ coronary artery disease,²⁶ tumor,^27,28 and other diseases. However, whether circRNAs show different expression profile in neuropathic pain states has not been reported.

To explore the expression of circRNAs in neuropathic pain condition, we developed the widely used CCI neuropathic pain model and detected the differentially expressed circRNAs between sham-operated and neuropathic pain rats with circRNA microarrays. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were carried out according to the microarray results and bioinformatics predictions and the circRNA–microRNA–mRNA interaction mode.

Methods

Animals

The experimental protocols were approved by the Experimental Animal Care and Use Committee of Zunyi Medical University. All experiments were conducted in accordance with the ethical guidelines of the International Association for the Study of Pain. All efforts were made to minimize animal suffering and to reduce animal use. Sprague Dawley rats (male, 7–8 weeks, 250–300 g) were housed under approved conditions with 12/12 hour light/dark cycles. Rats were free to eat and to drink water.

Induction of neuropathic pain

Rats underwent CCI surgeries to induce neuropathic pain on the left hind limb. Animals were randomly allocated to sham or CCI group. Neuropathic pain was induced with CCI as previously described by Bennett and Xie²⁹ with a 5-0 chromic gut suture.³⁰ In brief, after the intraperitoneal injection of sodium pentobarbital (35 mg/kg), the left sciatic nerve was loosely ligated with four sutures distant by 1–1.5 mm at the upstream of the nerve trifurcation of tibial, sural, and the common peroneal nerve. The muscle and skin were closed with sutures. Rats of the sham group underwent the same anesthesia and surgical procedures but the sciatic nerves were not ligated. Besides baseline tests, mechanical pain intensity was monitored at 2, 6, 10, and 14 days post-surgery.

Nociceptive behavioral testing

Mechanical hypersensitivity was determined by a male experimenter using an electronic von Frey plantar aesthesiometer (IITC, Wood Dale, IL, USA). Before tests, rats were given 10 minutes to habituate to the test environment. A rigid tip was applied against the mid-plantar surface of the left hind paw. The paw withdrawal threshold was automatically recorded by the device and the cutoff was set at 50 g. The rigid tip was presented perpendicularly to the plantar surface, and brisk withdrawal or paw flinching were considered as positive responses; the digital number presented on the monitor was recorded as the paw mechanical withdrawal threshold (MWT). Three successive stimuli were applied. MWT for individual animals was represented by the mean values (calculated from the three successive stimuli).

Ipsilateral dorsal horn samples preparation

After MWT tests on the 14th day post-surgery, rats were deeply anesthetized with isoflurane and were rapidly decapitated. The lumbar enlargement segments (L3–L5) of spinal cords were transversely sectioned and hemi-dissected along the midline. Only the dorsal half of the lumbar enlargement at the operated side of both sham-operated rats and CCI rats was collected and stored at −80°C.

RNA isolation, purification, and hybridization

One sample for each group (sham and CCI) was collected by pooling six ipsilateral L3–L5 dorsal spinal cords. RNA preparation and microarray hybridization were performed based on Arraystar’s standard protocols (Arraystar Inc, Rockville, MD, USA).³¹ The RNA amount of each sample was quantified using NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) and treated with Rnase R (Epicentre Inc, Madison, WI, USA) to remove linear RNAs and to enrich circRNAs. Then, the enriched circRNAs were amplified and transcribed into fluorescent cRNA using a random priming method (Arraystar Super RNA Labeling Kit; Arraystar Inc). The labeled cRNAs were hybridized onto the Arraystar Rat circRNA Array (8×15K, Arraystar Inc). At last, the arrays were scanned by an Agilent Scanner G2505C (Agilent Technologies, Santa Clara, CA, USA).

Microarray data analysis

Agilent Feature Extraction software (version 11.0; Agilent Technologies) was used to analyze acquired array images. Quantile normalization and subsequent data processing were performed using the R software Limma package (http://www.bioconductor.org/packages/release/bioc/html/limma.html). Differentially expressed circRNAs between two samples were identified through fold change filtering. CircRNAs having fold changes ≥2 were selected as of significantly differential expression. Hierarchical clustering was performed to show the distinguishable circRNAs expression pattern between sham and CCI samples.

Real-time polymerase chain reaction (PCR) validation

Quantitative real-time PCR (qPCR) was used to confirm the circRNA expression profiles obtained from the microarray data. Total RNA aforementioned was reversely transcribed into cDNA using a circRNA qRT-PCR kit (GeneSeed, Guangzhou, China) with TransScript RT enzyme (M-MLV) and reverse transcription primer according to the manufacturer’s protocol. The relative gene expression was determined using a CFX Connect Real-Time system (Bio-Rad Laboratories Inc, Hercules, CA, USA). All samples were normalized to the signal generated from GAPDH. Data are shown as fold change (2−ΔΔCt). Each sample was tested in triplicate. Primer sequences are listed in Table S1.

MicroRNA prediction

The circRNA–microRNA interaction was predicted with Arraystar’s home-made microRNA target prediction software based on the TargetScan³² and miRanda³³ prediction algorithm.

MicroRNA target genes prediction and bioinformatics analysis

To further investigate the functional roles of microRNA, putative targets of miRNAs were predicted by the TargetScan software. GO analysis was performed to explore the functional roles of target genes in terms of biological processes, cellular components, and molecular functions. Biological pathways defined by KEGG, Biocarta, and Reactome (http://www.genome.jp/kegg/) were identified by Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/).

CircRNA–microRNA interaction network

To further elucidate correlations between circRNAs and microRNA, potential circRNA–microRNA interaction analyses were conducted and maps were drawn with Cytoscape. The size of each node represents the number of putative microRNA functionally connected to each circRNA.

Statistical analysis

Results are reported as mean ± standard deviation. Statistically significant differences between the sham and CCI group were estimated by the repeated measures two-way analysis of variance using GraphPad Prism 7.0 (GraphPad Software Inc, La Jolla, CA, USA). P<0.05 was considered as being statistically significant.

Results

Development of mechanical hypersensitivity following CCI surgery

Compared with the sham group, rats in the CCI group showed significantly higher mechanical hypersensitivity on the ipsilateral hind paw at all of the four tested time points after CCI surgery (all P<0.001; Figure 1).

Figure 1 CCI-induced mechanical hypersensitivity.

Notes: Changes in MWT were assessed with the electronic von Frey filament on ipsilateral hind paws of CCI and sham-operated rats. n=6 in both groups. All data are represented as mean ± SD. Statistical analyses consisted of repeated measures two-way ANOVA tests. ***P<0.001 between sham and CCI groups.
Abbreviations: CCI, chronic constriction injury; SD, standard deviation; ANOVA, analysis of variance; MWT, mechanical withdrawal thresholds.

Overview of circRNA profiles

The Arraystar Rat circRNA Microarray detected 12,770 rat circRNAs (Supplementary materials). Hierarchical clustering and scatter plot visualization showed the circRNAs expression levels were different (Figure 2). Overall, 469 circRNAs were listed as differential circRNAs (fold change ≥2) between the CCI and sham group. Up to 363 circRNAs were significantly upregulated and 106 were downregulated in the CCI dorsal horn compared with that of sham-operated rats (Supplementary materials). The top 20 upregulated and downregulated circRNAs are listed in Table 1.

Figure 2 CircRNA expression between the CCI and sham-operated rats.

Notes: (A) Hierarchical clustering shows differential circRNA expression profiles between two groups. The upregulated circRNAs in the CCI group are shown in red, and the downregulated circRNAs are indicated in green. (B) Box plots show the distribution of circRNAs for the two compared samples. The distributions were nearly the same after normalization. (C) Scatter plots assess the circRNA expression differences between the two compared groups. The circRNAs above the top green line and below the bottom green line indicated differential circRNAs (fold change ≥2 or ≤−2, respectively).
Abbreviations: CCI, chronic constriction injury; circRNA, circular RNA.

Table 1 Top 20 upregulated and downregulated circRNAs between sham-operated and CCI rats (sham vs CCI)

CircRNA	Fold change	Best_transcript	Gene symbol	Raw intensity (CCI)	Raw intensity (sham)
Top 20 upregulated in CCI
rno_circRNA_013779	28.8	NM_001108050	Smek1	721	34
rno_circRNA_008008	11.3	XM_002728604	Nup133	1,517	139
rno_circRNA_003724	10.0	XM_222205	RGD1309762	1,257	128
rno_circRNA_016355	9.5	NM_021752	Birc2	1,451	158
rno_circRNA_008646	9.3	NM_001107689	Ash1l	3,225	333
rno_circRNA_011054	9.0	–	–	598	73
rno_circRNA_004560	8.7	NM_001100971	RGD1305110	818	104
mmu_circRNA_33592	8.6	NM_024486	Acvr1	978	122
rno_circRNA_007702	8.1	NM_001039713	Fto	4,271	489
rno_circRNA_001285	8.0	NM_001107604	Pdcd11	452	66
rno_circRNA_001148	7.9	NM_001107608	Hectd2	1,410	178
rno_circRNA_017183	7.9	NM_001191067	Tbc1d8	3,889	467
rno_circRNA_017427	7.7	NM_021687	Erbb4	1,301	174
rno_circRNA_009685	7.7	ENSRNOT00000071477	Mcu	1,158	155
rno_circRNA_007108	7.6	NM_001107398	Wdr33	582	85
rno_circRNA_008973	7.4	NM_017041	Ppp3ca	7,358	1,001
rno_circRNA_013293	7.4	NM_001107954	Reck	6,935	946
mmu_circRNA_000043	7.3	NM_001271502	Ctnnd2	995.5	145
rno_circRNA_014649	7.3	NM_001169103	Crim1	2,379	329
rno_circRNA_008686	7.1	NM_012780	Arnt	1,541	218
Top 20 downregulated in CCI
rno_circRNA_007512	−2.8	NM_001106133	Dym	138	382
rno_circRNA_007511	−2.6	NM_001106133	Dym	79	209
rno_circRNA_014170	−2.6	NM_001012191	Dtnb	272	703
rno_circRNA_016819	−2.6	ENSRNOT00000021758	Fert2	63	156
rno_circRNA_007514	−2.6	NM_001106133	Dym	90	228
rno_circRNA_007419	−2.5	NM_001107382	Zfp532	1,450	3,259
rno_circRNA_006343	−2.5	NM_031339	Parg	110	275
rno_circRNA_014174	−2.5	NM_001012191	Dtnb	272	666
rno_circRNA_009879	−2.5	NM_133300	Ddx39b	517	1,220
rno_circRNA_006094	−2.5	NM_001107312	Mtmr7	1,476	3,233
rno_circRNA_014172	−2.5	NM_001012191	Dtnb	245	589
rno_circRNA_014166	−2.5	NM_001012191	Dtnb	195	473
rno_circRNA_011941	−2.5	NM_030844	Ica1	188	458
rno_circRNA_002862	−2.5	NM_001100502	Gbe1	39	78
rno_circRNA_004376	−2.4	NM_138538	Dnm3	48	107
rno_circRNA_007562	−2.4	XM_225744	Setbp1	505	1,157
rno_circRNA_007421	−2.3	NM_001107382	Zfp532	872	1,939
rno_circRNA_013991	−2.3	NM_021767	Nrxn1	496.5	1,112
mmu_circRNA_45374	−2.3	XM_002730185	Kdm6a	76	177
rno_circRNA_001650	−2.3	NM_001191947	Dact3	327	764

Abbreviations: circRNAs, circular RNAs; CCI, chronic constriction injury.

qPCR validation

To validate the microarray data, qPCR was employed to detect circRNA expression. A total of eight circRNAs (four upregulated and four downregulated in the CCI group) were selected. Expression levels detected by the microarray and qPCR are presented in Figure 3. The expression trends detected by the two methods were consistent with each other, which demonstrated the high reliability of the microarray results.

Figure 3 Validation of eight differentially expressed circRNAs by qPCR.

Notes: qPCR results were calculated with 2−ΔΔCt method and normalized to GAPDH expression in the sham and CCI group. Bars represent the mean ± SD; positive fold change means overexpression in the CCI group. Each qPCR assay was performed three times.
Abbreviations: CCI, chronic constriction injury; SD, standard deviation; circRNA, circular RNA; qPCR, quantitative real-time polymerase chain reaction.

MicroRNA prediction and bioinformatics analyses

MicroRNA prediction was carried out with Arraystar’s homemade miRNA target prediction software based on miRanda and TargetScan. Target genes of the 505 predicted miRNAs of the top 100 differential circRNAs were predicted. We used miRDB (version 5) to forecast miRNAs’ target genes. Results of the GO enrichment analysis of the top 100 circRNAs with identified target genes are shown in Figure 4. GO analysis revealed thousands of target genes involved in the biological processes, such as cellular component organization (GO 0048856) etc, the cellular component, such as nucleoplasm (GO 0012505) etc, and molecular function, such as binding (GO 0005488) etc. KEGG analysis showed that there were 112 pathways related to the 100 circRNAs, including Hippo signaling pathway, MAPK signaling pathway, endocytosis and longevity regulating pathway etc (Supplementary materials). The top 20 high enrichment score pathways are listed in Figure 5.

Figure 4 GO enrichment analysis. GO enrichment corresponds to the top 100 circRNAs.

Abbreviations: circRNAs, circular RNAs; GO, gene ontology.

Figure 5 Top 20 KEGG pathways with the highest enrichment score.

Notes: Pathways correspond to the top 100 circRNAs. The vertical axis is the pathway category and the horizontal axis is the –LgP of the pathway; LgP is the logarithm of the P-value, and P<0.05 was considered significant.
Abbreviations: circRNAs, circular RNAs; KEGG, Kyoto Encyclopedia of Genes and Genomes.

CircRNA–microRNA competing endogenous RNAs (ceRNA) network

Potential connections between the top 10 circRNAs and microRNAs were predicted by using TargetScan and miRanda and the interactions were displayed with Cytoscape. As shown in Figure 6, seven of the top 10 circRNAs exhibited interactions with microRNAs. CircRNA_001059 and cir-cRNA_000167 were the two largest nodes in the network (with eight predicted microRNAs).

Figure 6 CircRNA–microRNA (ceRNA) network.

Notes: The circRNA–microRNA ceRNA network was drawn with the Cytoscape software. The size of each green node represents the functional connectivity of each circRNA. Rno_circRNA_008008 and rno_circRNA_013779 are the two largest nodes in the network.
Abbreviations: circRNAs, circular RNAs; ceRNA, competing endogenous RNAs.

Discussion

In the present study, we analyzed the circRNA expression in the dorsal horn of the spinal cord, which is an important location for pain regulation, synaptic plasticity, and neuropathic pain treatment.⁵ Comparison between sham-operated and CCI rats revealed expression change of 469 circRNAs. This study indicates that CCI of the peripheral nerve induces a lot of change of circRNAs in the rat spinal dorsal horn.

Despite the potential importance of circRNAs reported in several types of cancer,³⁴ heart diseases,^35–38 and other diseases, there is no report on the functional roles of circ-RNAs in the regulation of physical or pathological pain. In this study, in a rat neuropathic pain model, we found 363 circRNAs significantly upregulated and 106 circRNAs significantly downregulated in the spinal dorsal horn compared with normal ones. Three circRNAs (rno_circRNA_013779, rno_circRNA_013779, and rno_circRNA_013779) in CCI rats were upregulated more than 10 fold and the fold change of circRNA_013779 even reached 28. Interestingly, the number of downregulated circRNAs was less than that of the upregulated ones (363 vs 106); in addition, the magnitudes in the CCI group were much smaller. The top downregulated circRNA, rno_circRNA_007512, downregulated only 2.8 times in the CCI spinal dorsal horn. To evaluate the reliability of rat circRNA microarrays, we picked eight circRNAs of relatively big difference in the expression between two groups and with high expression levels in both groups. Expression patterns of these circRNAs were validated by qPCR. Consistent trends were observed between the microarrays and qPCR results (Figure 2).

Recent evidence revealed that circRNAs can function as microRNA sponges^22,39 and regulate parent gene expression to affect disease outcomes.^37,39 Despite the lack of literature reporting the exact functions of these differential circRNAs, we investigated the potential targets of the top 100 differential circRNAs (all of them were upregulated in the CCI group). In total, by using the miRanda software, 100 circRNAs were identified as the potential microRNA sponge of 505 micro-RNAs. GO analysis and KEGG pathway annotation were conducted to investigate the functions of related microR-NAs.³¹ GO enrichment analysis revealed that thousands of target genes were involved in the GO term: biological processes, cellular component, and molecular function, which indicates that these differential circRNAs may influence genes of varied biological function. KEGG analysis showed that there were 112 pathways related to the 100 circRNAs, the top three were reported to be related to neuropathic pain: the Hippo signaling pathway,⁴⁰ MAPK signaling pathway,^41–43 and endocytosis.^44,45

The circRNA–microRNA interaction network analysis was conducted for the top 10 circRNAs in this study. After filtering, we found that seven circRNAs may interact with microRNAs, all of them with more than three predicted target microRNAs. Two of them, circRNA_008008 and circRNA_013779, targeted eight microRNAs. CircRNAs were believed to negatively regulate miRNAs, and contribute substantially to the ceRNA network. It has been reported that ciRS-7, as a circular miR-7 inhibitor, harbors >60 conventional miR-7 binding sites, which is far more than any known linear sponge.⁴⁶ Sex determining region Y was another identified miRNA sponge and functioned as a miR-138 sponge.⁴⁷ According to our results, we hypothesized that some of these circRNAs may act as an inhibitor of miRNAs by acting as microRNA sponges and binding them. Our results implied that it is worthwhile to further investigate potential biological functions of these dysregulated circ-RNAs in the development of nerve constriction-induced neuropathic pain.

In the present study, although we pooled six dorsal horn samples in both groups, repeated microarray detections (ie, three repeats) and a set P threshold will be better to detect and decrease the number of differential circRNAs. In addition, comprehensive studies on circRNA expression at different time points and different nerve tissues (such as the peripheral sciatic nerve and dorsal root ganglion) in different neuropathic pain models are warranted. Literature has indicated that males and females respond differently in chronic pain conditions. Because the results of this study involve male rats,^48,49 we must realize that sex difference may influence the outcome of differential circRNAs. Functional studies of specific circRNAs in pain medicine are needed to clarify circRNAs roles in pain regulation.

Conclusion

Neuropathic pain rats displayed a circRNA expression change. Bioinformatics analyses showed that these circRNAs could be miRNA sponges and these microRNAs targeted mRNAs of various biological functions. CircRNA studies will be helpful to reveal molecular mechanisms of neuropathic pain.

Acknowledgments

The authors thank Dr Yi Zhang and Dr Wenjing Zhou for edits and revision of this manuscript. The study was supported by the National Natural Science Foundation of China (grant number 81660201 to Song Cao).

Disclosure

The authors report no conflicts of interest in this work.