Skip to main content
Download PDF

Systemic functional enrichment and ceRNA network identification following peripheral nerve injury

Molecular Brain volume 11, Article number: 73 (2018) Cite this article

Abstract

Peripheral nerve injury is a worldwide clinical issue that impacts patients’ quality of life and causes huge society and economic burden. Injured peripheral nerves are able to regenerate by themselves. However, for severe peripheral nerve injury, the regenerative abilities are very limited and the regenerative effects are very poor. A better understanding of the mechanisms following peripheral nerve injury will benefit its clinical treatment. In this study, we systematically explored the dynamic changes of mRNAs and long non-coding RNAs (lncRNAs) in the injured sciatic nerve segments after nerve crush, identified significantly involved Gene ontology (GO) terms and Kyoto Enrichment of Genes and Genomes (KEGG) pathways, and innovatively analyzed the correlation of differentially expressed mRNAs and lncRNAs. After the clustering of co-expressed mRNAs and lncRNAs, we performed functional analysis, selected GO term “negative regulation of cell proliferation”, and constructed a competing endogenous RNA (ceRNA) network of LIF and HMOX1 gene in this GO term. This study is the first to provide a systematic dissection of mRNA-microRNA (miRNA)-lncRNA ceRNA network following peripheral nerve injury and thus lays a foundation for further investigations of the regulating mechanisms of non-coding RNAs in peripheral nerve repair and regeneration.

Introduction

Peripheral nerve injury, resulting from a variety of different reasons such as mechanical compression, ischemia, penetrating injury, stretch injury, and cold injury, may seriously affect patients’ quality of life and cause huge society and economic burden [1]. It is reported that peripheral nerve injury affects about 1.5–2.8% of trauma patients and often leads to life-long morbidity and disability [2, 3]. In the United States, about 360,000 patients are suffering from upper extremity paralytic syndromes annually and more than $150 billion is used to treat peripheral nerve injury every year [4, 5].

Different from the central nervous system which has a very limited ability to regrowth, the peripheral nervous system obtains an intrinsic regenerative power [6]. In spite of this, the regenerative outcomes of injured peripheral nerves, especially peripheral nerves with severe defects and long nerve gaps, are generally poor and incomplete. Surgical nerve repair, including direct suturing and the transplantation of autologous nerve graft or tissue engineered nerve graft, improves the functional recovery of injured peripheral nerves [7]. However, the repairing effects of these therapeutical strategies are far from satisfactory. Gaining a better understanding of the cellular and molecular mechanisms underlying peripheral nerve injury and regeneration may contribute to the clinical treatment of peripheral nerve injury.

High-throughput screenings, such as microarrays and sequencing, are advanced large scale technologies for genome-wide analysis. Sequencing directly detects transcripts and has many advantages such as low background noise, large dynamic range, and high reproducibility [8, 9]. Besides the systematic identification of expressed profiles of known mRNAs, the application of sequencing also benefits the identifications of unannotated transcripts and non-coding RNAs, such as microRNAs (miRNAs), long non-coding RNAs (lncRNAs), and circular RNAs (circRNAs) [10,11,12,13]. In a previous study, by using RNA deep sequencing, we obtained the global transcriptome profiles of lesioned rat sciatic nerves at 0, 1, 4, 7, and 14 days after nerve crush. By using Ingenuity pathway analysis, we analyzed differentially expressed mRNAs and revealed key biological functions and canonical signaling pathways [14].

In the current study, with the joint use of Euclidean distance calculation, hierarchical clustering, principal component analysis, Gene ontology (GO), and Kyoto Enrichment of Genes and Genomes (KEGG), we further systematically determined the dynamic genetic changes following peripheral nerve injury. Besides a deeper investigation of differentially expressed mRNAs, we also identified differentially expressed lncRNAs, combined differentially expressed lncRNAs with differentially expressed mRNAs, and constructed correlated competing endogenous RNA (ceRNA) networks of LIF gene and HMOX1 gene based on miRWalk-validated miRNA-mRNA interaction and TargetScan-predicted miRNA-lncRNA interaction.

Materials and methods

RNA deep sequencing and data access

RNA sequencing was performed by using Illumina HiSeq™ 2000 and was described in the previous publication [14]. Briefly, rat sciatic nerve segments were collected from a total of 30 adult male Sprague-Dawley (SD) rats weighting 180–220 g at 0, 1, 4, 7, and 14 days after nerve injury. Total RNAs were extracted from sciatic nerve segments and purified. mRNAs and lncRNAs were fragmented into short pieces to synthesize cDNAs. cDNA fragments were purified, connected with adaptors, and used as templates for PCR amplification. Obtained raw reads subjected to quality control to collect clean reads by removing dirty reads with contain adaptors, high unknown bases, or low quality. Sequencing data were uploaded to NCBI database (accession number PRJNA394957; SRP113121).

Screening of significantly differentially expressed RNAs

The expression levels of both mRNAs and lncRNAs were calculated by using the Reads per kilobase transcriptome per million mapped reads (RPKM) method. The expression levels of mRNAs and lncRNAs at 1, 4, 7, and 14 days after rat sciatic nerve crush injury were compared to their expression levels at 0 day. RNAs with fold change> 10 and false discover rate (FDR) < 0.001 were screened and considered as significantly differentially expressed.

Bioinformatic analysis

Significantly differentially expressed RNAs at 1, 4, 7, and 14 days after nerve injury were analyzed by using the Venny 2.1.0 software (http://bioinfogp.cnb.csic.es/tools/venny/index.html) [15, 16] to visualize the intersection of differentially expressed RNAs at each time point. Euclidean distance calculation and hierarchical clustering were performed by using the HeatMapImage GenePattern module to illustrate the temporal expression patterns of differentially expressed RNAs. Principal component analysis was performed by using the Population Principal Component Analysis software (Harvard Medical School, MA, USA) to display the similarity of differentially expressed RNAs at different time point. Database for Annotation, Visualization, and Integrated Discovery (DAVID) bioinformatic resource (https://david.ncifcrf.gov/) was used to identify enriched GO terms and KEGG pathways.

Significantly differentially expressed mRNAs and lncRNAs were subjected to the calculation of Pearson correlation coefficient. mRNAs and lncRNAs with Pearson correlation coefficient index> 0.9 and adjusted p-value< 0.1 were considered as co-expressed. K-means clustering was conducted and GO and KEGG analysis was performed to enrich GO terms and KEGG pathways with p-value< 0.05 in each cluster. The binding relationships of mRNAs and miRNAs were analyzed and validated by using the miRWalk software (http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/gopub.html). GO terms with more than one validated gene were screened and GO term “negative regulation of cell proliferation” was selected for the construction of ceRNA network. The binding relationships of miRNAs and lncRNAs were predicted by using TargetScan software (http://www.targetscan.org/vert_72/). The interactions of mRNAs LIF, HMOX1, validated miRNAs miR-494-3p, let-7e-5p, let-7a-5p, let-7d-5p, and predicted lncRNAs were analyzed and corresponding ceRNA network was built.

Animal surgery

To validate outcomes from bioinformatic analysis, we obtained a total of 45 adult male SD rats (180–220 g) from the Experimental Animal Center of Nantong University and performed sciatic nerve crush injury as previously described [14]. Briefly, after anaesthetization, rat sciatic nerve at 10 mm above the bifurcation into the tibial and common fibular nerves was crushed with a forceps at a force of 54 N for 30 s with 10 s for each time. Rats were sacrificed by cervical dislocation at 1, 4, 7, and 14 days after surgery. Rats underwent sham-surgery were used as 0 day control. All the experimental procedures involving animals were conducted in accordance with Institutional Animal Care guidelines of Nantong University and approved ethically by the Administration Committee of Experimental Animals, Jiangsu, China.

Quantitative RT-PCR

Sciatic nerve segments at the crush site (different samples from those for RNA deep sequencing) were collected to extract total RNAs by using Trizol Reagent (Life Technologies, Carlsbed, CA, USA). Isolated RNAs were reverse transcribed to cDNAs by using Prime-Script Reagent Kit (TaKaRa, Dalian, Liaoning, China) or TaqMan MicroRNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA, USA). cDNAs were amplified with QuantiNova SYBR Green PCR Kit (Qiagen, Valencia, CA, USA) on an Applied Biosystems Stepone real-time PCR System (Applied Biosystems, Foster City, CA, USA). Relative abundances of mRNAs, miRNAs, and lncRNAs were determined by using the ΔΔCt method. GAPDH was used as the reference gene for the quantifications of mRNAs and lncRNAs and U6 was used as the reference gene the quantifications of miRNAs. Primers used were listed in Additional file 1: Table S1.

Statistical analysis

Statistical analyses were performed by using GraphPad Prism 6.0 (GraphPad Software, Inc., La Jolla, CA, USA). One-way analysis of variance (ANOVA) and Dunnett’s multiple comparisons test were used for comparison between groups. Numerical data were presented as mean ± SEM and a p-value< 0.05 was considered as statistically different.

Results

Identification of significantly differentially expressed RNAs after sciatic nerve injury

Our previously obtained sequencing data (SRP113121) discovered a total of 38,967 RNAs (including mRNAs and lncRNAs) in rat sciatic nerve segments with 35,728, 38,024, 36,847, 37,513, and 36,403 RNAs at 0, 1, 4, 7, and 14 days, respectively [14]. After the comparison of the expressions of RNAs in injured rats with those at 0 day, a total of 22,498 RNAs were found to be differentially expressed during the time course (fold change> 2 or < − 2, FDR < 0.001). The vast numerical RNAs increased the difficulty of the discovery of critical information. Therefore, we filtered sequencing data, selected RNAs with fold change> 10 or < − 10 and FDR < 0.001 as compared with 0 day control, and designated these RNAs as significantly differentially expressed RNAs. A full list of all significantly differentially expressed mRNAs and lncRNAs were shown in Additional file 2: Table S2.

It was demonstrated that compared with 0 day control, at 1 day after nerve injury, 575 mRNAs and 295 lncRNAs were up-regulated while 70 mRNAs and 17 lncRNAs were down-regulated. Slightly smaller numbers of RNAs were significantly differentially expressed at 4, 7, and 14 days after nerve injury. At 14 day after nerve injury, only 354 RNAs (240 mRNAs and 114 lncRNAs) were up-regulated and 90 RNAs (80mRNAs and 10 lncRNAs) were down-regulated (Fig. 1a). Significantly differentially expressed RNAs during the time course were further analyzed and illustrated by the Venn diagram. A total of 166 RNAs were commonly up-regulated or down-regulated at all time points while many RNAs were only significantly differentially expressed at one single time point (Fig. 1b). Differentially expressed RNAs were then subjected to cluster analysis to determine the similarity of RNA expression profiles during the time course. Euclidean distance, hierarchical clustering, and principal component analysis demonstrated that the RNA expression patterns in 0 and 1 day were quite different while the RNA expression patterns in 4, 7, and 14 day after nerve injury showed certain consistency (Fig. 1c & d).

Fig. 1
figure 1

Overview of significantly differentially expressed RNAs in the sciatic nerve segments after injury. a The bar graph of significantly differentially expressed mRNAs and lncRNAs at 1, 4, 7, and 14 days after sciatic nerve injury. RNAs with fold change> 10 or < − 10 and FDR < 0.001 as compared with 0 day control were considered as significantly differentially expressed. DEG, differentially expressed genes. Red color indicated up-regulated RNAs while green color indicated down-regulated RNAs. b The Venn diagram of significantly differentially expressed mRNAs and lncRNAs at 1, 4, 7, and 14 days after sciatic nerve injury. c Euclidean distance and hierarchical clustering of RNAs at 0, 1, 4, 7, and 14 days after sciatic nerve injury. d Principal component analysis of RNAs at 0, 1, 4, 7, and 14 days after sciatic nerve injury

Full size image

Functional enrichment analysis of differentially expressed mRNAs

Significantly differentially expressed mRNAs were then categorized to GO terms and KEGG pathways to discover possible biological functions and pathways involved in peripheral nerve injury. A full list of GO biological process, cellular component, and molecular function terms with a p-value< 0.05 at 1, 4, 7, and 14 days after nerve injury was provided in Additional file 3: Table S3. GO terms with a p-value< 0.05 and FDR < 0.05 at least one time point were enriched, listed, ranked with p-value and FDR from day 1 to day 14, and presented in color gradation (Fig. 2). It was observed that mass GO terms were dramatically involved at 1 day after nerve injury while less GO terms were significantly enriched at later time points. For GO cellular component terms, extracellular space, cell surface, and external side of plasma membrane were the top three enriched terms. For GO molecular function terms, chemokine activity, cytokine activity, and carbohydrate binding were the top three enriched terms. A larger number of GO biological process terms were involved. Notably, many of these enriched GO biological process terms were associated with immune and inflammatory response (neutrophil chemotaxis, inflammatory response, immune response, cellular response to interleukin-1, lymphocyte chemotaxis, cellular response to interferon-gamma, defense response to Gram-positive bacterium, defense response to bacterium, positive regulation of inflammatory response, cellular response to tumor necrosis factor, positive regulation of neutrophil chemotaxis, acute inflammatory response, monocyte chemotaxis, innate immune response, leukocyte cell-cell adhesion, leukocyte chemotaxis, response to bacterium, and positive regulation of T cell proliferation). Many cell signaling pathway-related GO biological process terms, for example, chemokine-mediated signaling pathway, cytokine-mediated signaling pathway, positive regulation of ERK1 and ERK2 cascade, were also enriched.

Fig. 2
figure 2

Enriched GO terms of significantly differentially expressed mRNAs in the sciatic nerve segments after injury. GO cellular component, molecular function, and biological process terms with p-value< 0.05 and FDR < 0.05 in 1, 4, 7, or 14 days after sciatic nerve injury were screened. The p-values of GO terms were listed. GO terms with low p-values were labeled in red color while GO terms with high p-values were labeled in green color

Full size image

To further investigate the involvement of cell signaling pathways in peripheral nerve injury, KEGG pathway analysis was conducted and critical signaling pathways in peripheral nerve injury were identified (Fig. 3). A total of seven signaling pathways, including chemokine signaling pathway, cytokine-cytokine receptor interaction, hematopoietic cell lineage, malaria, neuroactive ligand-receptor interaction, osteoclast differentiation, and rheumatoid arthritis, were significantly activated at all time points after nerve injury. And signaling pathways chemokine signaling pathway, cytokine-cytokine receptor interaction, and neuroactive ligand-receptor interaction obtained large gene numbers and high Rich Factor, implying their essential roles in peripheral nerve injury.

Fig. 3
figure 3

Enriched KEGG pathways of significantly differentially expressed mRNAs in the sciatic nerve segments after injury. KEGG pathways with p-value< 0.05 in 1, 4, 7, or 14 days after sciatic nerve injury were screened and listed. The number of significantly differentially expressed genes in each KEGG pathway and relevant Rich Factor were shown

Full size image

Co-expression of significantly differentially expressed mRNAs and lncRNAs

Besides the functional analysis of mRNAs, we further determined the correlation of significantly differentially expressed mRNAs and lncRNAs, aiming to characterize mRNA-miRNA-lncRNA ceRNA networks in peripheral nerve injury (Fig. 4a). Person correlation analysis showed that many significantly differentially expressed mRNAs and lncRNAs were correlated with each other (Fig. 4b). Co-expressed mRNAs and lncRNAs were then subjected to K-means clustering to group RNAs with different temporal expression patterns. Two clusters of RNAs (boxed by different colors) were categorized and functional enriched to GO terms and KEGG pathways to discover critical biological activities of co-expressed mRNAs and lncRNAs (Fig. 4c). By setting a cutoff of p-value< 0.05, we identified a total of 317 activated GO terms and KEGG pathways in cluster 1 and a total of 130 activated GO terms and KEGG pathways in cluster 2 (Additional file 4: Table S4). We further studied involved mRNAs in these GO terms and KEGG pathways and determined the validated binding relationships between mRNAs and miRNAs by using the miRWalk software. A total of 54 GO terms in cluster 1 and 10 GO terms in cluster 2 were found to have mRNAs with validated miRNA binding relationships (Additional file 5: Table S5). Enriched GO terms with more than one validated mRNA were screened and listed (Fig. 4c).

Fig. 4
figure 4

Coexpression analysis of significantly differentially expressed mRNAs and lncRNAs in the sciatic nerve segments after injury. a Schematic depiction of coexpression analysis of significantly differentially expressed mRNAs and lncRNAs and the construction of mRNA-miRNA-lncRNA ceRNA network. b Person correlation network of significantly differentially expressed mRNAs and lncRNAs. mRNAs were labeled in cyan color while lncRNAs were labeled in red color. c K-means clustering of co-expressed mRNAs and lncRNAs. Enriched GO terms and involved validated genes were listed

Full size image

Construction and examination of LIF and HMOX1-associated ceRNA network

Enriched GO term “negative regulation of cell proliferation” contained three validated genes leukemia inhibitory factor (LIF), heme oxygenase (decycling) 1 (HMOX1), and peripheral myelin protein 22 (PMP22). Therefore validated genes in this GO term were further investigated. Outcomes from miRWalk database suggested that in GO term “negative regulation of cell proliferation”, LIF interacted with miR-494-3p, HMOX1 interacted with miR-494-3p, let-7e-5p, let-7a-5p, and let-7d-5p, and PMP22 interacted with miR-9a-5p and miR-29a-3p. Since both LIF and HMOX1 interacted with miR-494-3p, a network containing LIF, HMOX1, and their validated binding miRNAs was built (Fig. 5a). Moreover, we predicted lncRNAs that bound to miRNAs in the network by TargetScan, jointly analyzed mRNA-miRNA-lncRNA interactions, and constructed a LIF and HMOX1-associated ceRNA network (Fig. 5a). RNAs in the ceRNA network were further investigated. Sequencing outcomes suggested that both LIF and HMOX1 were up-regulated after nerve injury (Fig. 5b). The abundances of lncRNAs were also determined and demonstrated in heatmaps (Fig. 5c).

Fig. 5
figure 5

ceRNA network of LIF and HMOX1. a Interactions of RNAs in LIF and HMOX1-associated ceRNA network. b Expression levels of mRNAs LIF and HMOX1 from RNA sequencing outcome. c Heatmap and hierarchical clustering of lncRNAs in LIF and HMOX1-associated ceRNA network. Red color indicated up-regulation while green color indicated down-regulation

Full size image

The temporal expression levels of RNAs in the ceRNA network were further determined by quantitative RT-PCR. RT-PCR results showed that consistent with sequencing outcomes, the expression levels of LIF were increased after nerve injury (Fig. 6a). In contrast, the expression levels of miR-494-3p were decreased after nerve injury (Fig. 6b). The negative correlation of the expression patterns of LIF and miR-494-3p further demonstrated that LIF might be the mRNA target of miR-494-3p in the sciatic nerve segments after peripheral nerve injury. The expression levels of lncRNAs XLOC_083376, XLOC_150234, XLOC_138166, XLOC_027762, XLOC_134667, XLOC_056487, XLOC_146842, XLOC_174535, XLOC_174536, XLOC_065278, and XLOC_083385 were also determined. PCR results showed that XLOC_083376, XLOC_150234, XLOC_138166, XLOC_134667, XLOC_146842, XLOC_174535, XLOC_174536, XLOC_065278, and XLOC_083385 were mainly up-regulated after nerve injury, XLOC_056487, and XLOC_065278 were down-regulated, while XLOC_027762 was first down-regulated and then up-regulated (Fig. 6c–m).

Fig. 6
figure 6

PCR validation of RNAs in the LIF-associated ceRNA network. The expression levels of (a) LIF, (b) miR-494-3p, (c) XLOC_083376, (d) XLOC_150234, (e) XLOC_138166, (f) XLOC_027762, (g) XLOC_134667, (h) XLOC_056487, (i) XLOC_146842, (j) XLOC_174535, (k) XLOC_174536, (l) XLOC_065278, and (m) XLOC_083385 at 0, 1, 4, 7, and 14 days after sciatic nerve injury were determined and normalized to reference RNA GAPDH. Numerical data were summarized from three independent experiments. *p-value< 0.05

Full size image

Similarly, the abundances of RNAs in the HMOX1-let-7 ceRNA network were also determined. Following sciatic nerve injury, the expression levels of HMOX1 and let-7 (let-7e, let-7a, and let-7d) were increased and decreased, respectively (Fig. 7a–d). Together with observation in Fig. 6b, RT-PCR outcomes demonstrated that temporal expression changes of HMOX1 were inversely associated with those of miR-494-3p, let-7e, let-7a, and let-7d. Determination of lncRNA expression patterns showed that XLOC_080446, XLOC_174287, XLOC_000412, XLOC_130885, and XLOC_150864 were up-regulated, XLOC_172336, XLOC_146853, and XLOC_133837 were down-regulated, while XLOC_150864 and XLOC_108671 were first down-regulated and then up-regulated expressed.

Fig. 7
figure 7

PCR validation of RNAs in the HMOX1-associated ceRNA network. The expression levels of (a) HMOX1, (b) let-7e-5p, (c) let-7a-5p, (d) let-7d-5p, (e) XLOC_080446, (f) XLOC_172336, (g) XLOC_174287, (h) XLOC_000412, (i) XLOC_130885, (j) XLOC_146853, (k) XLOC_150864, (l) XLOC_133837, and (m) XLOC_108671 at 0, 1, 4, 7, and 14 days after sciatic nerve injury were determined and normalized to reference RNA GAPDH. Numerical data were summarized from three independent experiments. *p-value< 0.05

Full size image

Discussion

In the current study, we used RNA deep sequencing and bioinformatic analysis to investigate molecular changes following peripheral nerve injury. Previously, by comparing the expressions of RNAs at different time points to their expressions at 0 day and by setting a threshold of fold change> 2 or < − 2 and FDR < 0.001, we identified 13,721, 14,321, 14,745, and 6979 differentially expressed RNAs at 1, 4, 7, and 14 days after nerve crush [14]. To screen out RNAs with extreme changes, here, we increased the threshold to fold change> 10 or < − 10 and FDR < 0.001 and discovered 957, 886, 590, and 444 significantly differentially expressed RNAs at 1, 4, 7, and 14 days after nerve crush, respectively. Moreover, we separately counted the numbers of differentially expressed mRNAs and lncRNAs and found that about 1/3 of differentially expressed RNAs were lncRNAs. The large amount of differentially expressed lncRNAs may further affect tons of mRNAs since lncRNAs could modulate the expressions of many mRNAs in multiple levels, including transcriptional regulation, post-transcriptional regulation, and epigenetic modification [17,18,19].

Euclidean distance, hierarchical clustering, and principal component analysis outcomes demonstrated that RNA expressions in the uninjured 0 day group were distinct from RNA expressions in the lesioned sciatic nerves. A comparison of RNA expressions at different time points after sciatic nerve injury showed that the expression profiles of RNAs in 1 day were also obviously different from those at later time points. GO annotation showed that a large number of GO cellular component, molecular function, and biological process terms were significantly enriched at 1 day after nerve crush while relatively smaller numbers of GO terms were enriched at later time points. Some inflammatory and immune response-related GO biological process terms (e.g., neutrophil chemotaxis, inflammatory response, and immune response) were kept activated at all time points. This was consistent with our previous observations of the importance of inflammatory and immune response after nerve injury [14, 20]. Enriched signaling pathways were also examined by KEGG pathway analysis. Our previous microarray analysis of the distal sciatic nerve segments showed that cytokine-cytokine receptor interaction and neuroactive ligand-receptor interaction were critical signaling pathways during Wallerian degeneration [21,22,23,24,25]. Here, we found that in the entire lesioned nerve segments, these two signaling pathways were also significantly involved. Some other signaling pathways, for example, chemokine signaling pathway, were demonstrated to be activated in the entire lesioned nerve segments but not the distal nerve segments. It implied that these signaling pathways might be important for nerve regrowth from the proximal nerve segments.

Besides these bioinformatic analyses, in the current study, we also jointly analyzed differentially expressed mRNAs and lncRNAs, identified the correlation of differentially expressed mRNAs and lncRNAs, clustered co-expressed mRNAs and lncRNAs, and performed functional analysis of co-expressed mRNAs and lncRNAs. Enriched GO terms and KEGG pathways in each cluster were identified and functional terms with validated mRNAs were further selected for the construction of ceRNA networks. There existed eight enriched GO terms (negative regulation of cell proliferation, cytosol, positive regulation of angiogenesis, response to hypoxia, response to nicotine, extracellular space, plasma membrane, and protein homodimerization activity) with more than one validated mRNAs in cluster 1 and one enriched GO term (neuron projection) with more than one validated mRNAs in cluster 2. We then selected GO term “negative regulation of cell proliferation”, identified bound miRNAs of validated genes in the GO terms by using miRWalk database, predicted interacted lncRNAs by using TargetScan software, and constructed LIF and HMOX1-associated ceRNA network.

Emerging studies have shown that miRNAs are key regulators in many physiological and pathological processes [26,27,28]. It has been demonstrated that many miRNAs were differentially expressed after peripheral nerve injury [29, 30]. Dysregulated miRNAs regulate Schwann cell proliferation, migration, and myelination and affect peripheral nerve regeneration [31,32,33,34]. Notably, the biological effects of miRNAs can be modulated by mRNAs, transcribed pseduogenes, lncRNAs, and circRNAs through the competitive binding to the miRNA response elements [35,36,37,38].

As far as we know, till now, there is no study about lncRNA-associated ceRNAs in peripheral nerve repair and regeneration. Here, we used computational methods to discover correlations between mRNAs and lncRNAs and used miRWalk and TargetScan algorithm to explore mRNA-miRNA-lncRNA interactions. The temporal expression levels of mRNAs and lncRNAs in the constructed LIF and HMOX1-associated ceRNA network were also determined and demonstrated in line charts and heatmaps. Furthermore, the expression levels of mRNAs, miRNAs, and lncRNAs in the ceRNA network were validated by RT-PCR. Sequencing and RT-PCR results showed that both the expressions of LIF and HMOX1 were up-regulated after nerve injury. In contrast, the expressions of miR-494-3p, let-7e-5p, let-7a-5p, and let-7d-5p, validated bound miRNAs of LIF and HMOX1, were down-regulated. This, from the aspect of RNA expression, further supported that miR-494-3p, let-7e-5p, let-7a-5p, and let-7d-5p might regulate LIF and HMOX1 after peripheral nerve injury. Moreover, RT-PCR results showed that lncRNAs XLOC_174535, XLOC_174536, and XLOC_000412 were up-regulated after nerve injury, this was consistent with the expressions of LIF and HMOX1 and opposite with the expressions of miR-494-3p, let-7e-5p, let-7a-5p, and let-7d-5p. The correlations between the expressions of mRNAs and these lncRNAs were further calculated by linear regression (Fig. 8). The R2 of LIF expression and XLOC_174535 expression was 0.996 and the R2 of LIF expression and XLOC_174536 expression was 0.708 (Fig. 8a). The high R2 value suggested that LIF is positively correlated with lncRNAs XLOC_174535 and XLOC_174536, indicating that in LIF ceRNA network, lncRNAs XLOC_174535 and/or XLOC_174536 might sponge miR-494-3p and regulate LIF expression. Similarly, linear regression calculation showed that the R2 of HMOX1 expression and XLOC_174535 expression was 0.851 while the R2 of HMOX1 expression and XLOC_174535 expression was a little bit lower (0.581) (Fig. 8b). Therefore, in the HMOX1-miR-494-3p-lncRNA ceRNA network, lncRNA XLOC_174535 might sponge miR-494-3p and regulate HMOX1 expression. In HMOX1-let-7e-5p/let-7a-5p/let-7d-5p-lncRNA ceRNA network, the R2 of HMOX1 and XLOC_000412 was 0.687 (Fig. 8b), telling that lncRNA XLOC_000412 might sponge let-7e-5p/let-7a-5p/let-7d-5p and regulate HMOX1 expression. Luciferase assay could be performed to further examine the relationships of RNAs in the ceRNA network.

Fig. 8
figure 8

Correlation of mRNA and lncRNA expressions. a Expression correlation of LIF and lncRNAs XLOC_174535 and XLOC_174536. b Expression correlation of HMOX1 and lncRNAs XLOC_174535, XLOC_174536, and XLOC_000412. R2 was calculated by linear regression

Full size image

On the other hand, it is worth noting that our validation outcomes showed that the expression profiles of some lncRNAs were conflicting with sequencing outcomes. Many factors, including adequate replication, the quality of RNA, different efficiencies of reverse transcriptases, varied priming methods, and data normalization differences, can affect quantitative results and may lead to inconsistent results between outcomes from high-throughput analysis and RT-PCR validation [39]. This inconsistency was also observed in other studies [40,41,42]. Morey et al. calculated the direction of change in expression (up-regulation or down-regulation) by microarray and PCR and showed that the direction of change in expression was in agreement for 72.9% of samples [39]. Here, we found that in all tested 20 lncRNAs, 5 lncRNAs (XLOC_056487, XLOC_065278, XLOC_172336, XLOC_146853, and XLOC_133837) showed conflicting directions (down-regulated instead of up-regulated after nerve injury). The ratio of the direction of change in expression (25% disagreement) in our current study was similar as Morey’s calculation. In addition, we examined the sequencing RPKM reads of tested mRNAs and lncRNAs (Additional file 6: Table S6) and found that compared with the PRKM reads of HMOX1 and LIF, the PRKM reads of many lncRNAs were much lower. Genes with low absolute expression levels normally had larger changes of having inconsistent validation results [39]. Therefore, it is possible that the accuracy of lncRNAs with lower RPKM reads may be influenced by their lower expression levels.

In summary, in the current study, we studied the temporal changes of mRNAs and lncRNAs in the lesioned sciatic nerve segments at 0, 1, 4, 7, and 14 days after nerve crush and detected enriched GO terms and KEGG pathways. Moreover, for the first time, we elucidated the correlations of differentially expressed mRNAs and lncRNAs and delineated functional landscapes of mRNA-miRNA-lncRNA ceRNA network following peripheral nerve injury. Our current study expanded our knowledge about the molecular basis of peripheral nerve injury, provided insights of the potential regulations of non-coding RNAs, and offered promising prospects of the clinical treatment of peripheral nerve injury.

Abbreviations

ANOVA:

One-way analysis of variance

ceRNA:

Competing endogenous RNA

circRNA:

Circular RNA

DAVID:

Database for Annotation, Visualization, and Integrated Discovery

FDR:

False discover rate

GO:

Gene ontology

HMOX1:

Heme oxygenase (decycling) 1

KEGG:

Kyoto Enrichment of Genes and Genomes

LIF:

Leukemia inhibitory factor

lncRNA:

Long non-coding RNA

miRNA:

MicroRNA

PMP22:

Peripheral myelin protein 22

RPKM:

Reads per kilobase transcriptome per million mapped reads

SD:

Sprague-Dawley

References

  1. Gu X, Ding F, Yang Y, Liu J. Construction of tissue engineered nerve grafts and their application in peripheral nerve regeneration. Prog Neurobiol. 2011;93(2):204–30. https://doi.org/10.1016/j.pneurobio.2010.11.002.

    Article  CAS  PubMed  Google Scholar 

  2. Noble J, Munro CA, Prasad VS, Midha R. Analysis of upper and lower extremity peripheral nerve injuries in a population of patients with multiple injuries. J Trauma. 1998;45(1):116–22.

    Article  CAS  Google Scholar 

  3. Taylor CA, Braza D, Rice JB, Dillingham T. The incidence of peripheral nerve injury in extremity trauma. Am J Phys Med Rehabil. 2008;87(5):381–5. https://doi.org/10.1097/PHM.0b013e31815e6370.

    Article  PubMed  Google Scholar 

  4. Tian L, Prabhakaran MP, Ramakrishna S. Strategies for regeneration of components of nervous system: scaffolds, cells and biomolecules. Regen Biomater. 2015;2(1):31–45. https://doi.org/10.1093/rb/rbu017.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Barton N. Upper extremity disorders: frequency, impact and cost, J. Kelsey, A. Praemer, L. Nelson, A. Felberg, D. Rice. Churchill Livingstone, New York (1997), Price £15.95, paperback, ISBN: 0-443-07912-9. J Hand Surg Br Eur. 1998;23(2):285.

    Article  Google Scholar 

  6. Chen ZL, Yu WM, Strickland S. Peripheral regeneration. Annu Rev Neurosci. 2007;30:209–33. https://doi.org/10.1146/annurev.neuro.30.051606.094337.

    Article  CAS  PubMed  Google Scholar 

  7. Li R, Liu Z, Pan Y, Chen L, Zhang Z, Lu L. Peripheral nerve injuries treatment: a systematic review. Cell Biochem Biophys. 2014;68(3):449–54. https://doi.org/10.1007/s12013-013-9742-1.

    Article  CAS  PubMed  Google Scholar 

  8. Nagalakshmi U, Waern K, Snyder M. RNA-Seq: a method for comprehensive transcriptome analysis. Curr Protoc Mol Biol. 2010;Chapter 4(Unit 4):4.11.11–14.11.13.

    Google Scholar 

  9. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8. https://doi.org/10.1038/nmeth.1226.

    Article  CAS  Google Scholar 

  10. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5. https://doi.org/10.1038/nbt.1621.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Xu T, Wu J, Han P, Zhao Z, Song X. Circular RNA expression profiles and features in human tissues: a study using RNA-seq data. BMC Genomics. 2017;18(Suppl 6):680. https://doi.org/10.1186/s12864-017-4029-3.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Fan J, Zhou Q, Li Y, Song X, Hu J, Qin Z, Tang J, Tao T. Profiling of long non-coding RNAs and mRNAs by RNA-sequencing in the hippocampi of adult mice following propofol sedation. Front Mol Neurosci. 2018;11:91. https://doi.org/10.3389/fnmol.2018.00091.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Guelfi G, Cochetti G, Stefanetti V, Zampini D, Diverio S, Boni A, Mearini E. Next Generation Sequencing of urine exfoliated cells: an approach of prostate cancer microRNAs research. Sci Rep. 2018;8(1):7111. https://doi.org/10.1038/s41598-018-24236-y.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Yi S, Zhang H, Gong L, Wu J, Zha G, Zhou S, Gu X, Yu B. Deep sequencing and bioinformatic analysis of lesioned sciatic nerves after crush injury. PLoS One. 2015;10(12):e0143491. https://doi.org/10.1371/journal.pone.0143491.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Chen H, Boutros PC. VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics. 2011;12:35. https://doi.org/10.1186/1471-2105-12-35.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Oliveros JC. VENNY. An interactive tool for comparing lists with Venn Diagrams; 2007.

    Google Scholar 

  17. Guttman M, Rinn JL. Modular regulatory principles of large non-coding RNAs. Nature. 2012;482(7385):339–46. https://doi.org/10.1038/nature10887.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Wang X, Song X, Glass CK, Rosenfeld MG. The long arm of long noncoding RNAs: roles as sensors regulating gene transcriptional programs. Cold Spring Harb Perspect Biol. 2011;3(1):a003756. https://doi.org/10.1101/cshperspect.a003756.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kurokawa R, Rosenfeld MG, Glass CK. Transcriptional regulation through noncoding RNAs and epigenetic modifications. RNA Biol. 2009;6(3):233–6.

    Article  CAS  Google Scholar 

  20. Xing L, Cheng Q, Zha G, Yi S. Transcriptional profiling at high temporal resolution reveals robust immune/inflammatory responses during rat sciatic nerve recovery. Mediat Inflamm. 2017;2017:3827841. https://doi.org/10.1155/2017/3827841.

    Article  CAS  Google Scholar 

  21. Yi S, Tang X, Yu J, Liu J, Ding F, Gu X. Microarray and qPCR analyses of Wallerian degeneration in rat sciatic nerves. Front Cell Neurosci. 2017;11:22. https://doi.org/10.3389/fncel.2017.00022.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Cheng Q, Wang YX, Yu J, Yi S. Critical signaling pathways during Wallerian degeneration of peripheral nerve. Neural Regen Res. 2017;12(6):995–1002. https://doi.org/10.4103/1673-5374.208596.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Li M, Zhang P, Guo W, Li H, Gu X, Yao D. Protein expression profiling during wallerian degeneration after rat sciatic nerve injury. Muscle Nerve. 2014;50(1):73–8. https://doi.org/10.1002/mus.24082.

    Article  CAS  PubMed  Google Scholar 

  24. Yao D, Li M, Shen D, Ding F, Lu S, Zhao Q, Gu X. Expression changes and bioinformatic analysis of Wallerian degeneration after sciatic nerve injury in rat. Neurosci Bull. 2013;29(3):321–32. https://doi.org/10.1007/s12264-013-1340-0.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Li M, Guo W, Zhang P, Li H, Gu X, Yao D. Signal flow and pathways in response to early Wallerian degeneration after rat sciatic nerve injury. Neurosci Lett. 2013;536:56–63. https://doi.org/10.1016/j.neulet.2013.01.008.

    Article  CAS  PubMed  Google Scholar 

  26. Fitzpatrick JM, Anderson RC, McDermott KW. MicroRNA: key regulators of oligodendrocyte development and pathobiology. Int J Biochem Cell Biol. 2015;65:134–8. https://doi.org/10.1016/j.biocel.2015.05.021.

    Article  CAS  PubMed  Google Scholar 

  27. Herkenhoff ME, Oliveira AC, Nachtigall PG, Costa JM, Campos VF, Hilsdorf AWS, Pinhal D. Fishing into the MicroRNA transcriptome. Front Genet. 2018;9:88. https://doi.org/10.3389/fgene.2018.00088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Geng L, Sun B, Gao B, Wang Z, Quan C, Wei F, Fang XD. MicroRNA-103 promotes colorectal cancer by targeting tumor suppressor DICER and PTEN. Int J Mol Sci. 2014;15(5):8458–72.

    Article  CAS  Google Scholar 

  29. Li S, Yu B, Wang Y, Yao D, Zhang Z, Gu X. Identification and functional annotation of novel microRNAs in the proximal sciatic nerve after sciatic nerve transection. Sci China Life Sci. 2011;54(9):806–12. https://doi.org/10.1007/s11427-011-4213-7.

    Article  CAS  PubMed  Google Scholar 

  30. Yu B, Zhou S, Wang Y, Ding G, Ding F, Gu X. Profile of microRNAs following rat sciatic nerve injury by deep sequencing: implication for mechanisms of nerve regeneration. PLoS One. 2011;6(9):e24612. https://doi.org/10.1371/journal.pone.0024612.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Li S, Zhang R, Yuan Y, Yi S, Chen Q, Gong L, Liu J, Ding F, Cao Z, Gu X. MiR-340 regulates fibrinolysis and axon regrowth following sciatic nerve injury. Mol Neurobiol. 2017;54(6):4379–89. https://doi.org/10.1007/s12035-016-9965-4.

    Article  CAS  PubMed  Google Scholar 

  32. Yi S, Yuan Y, Chen Q, Wang X, Gong L, Liu J, Gu X, Li S. Regulation of Schwann cell proliferation and migration by miR-1 targeting brain-derived neurotrophic factor after peripheral nerve injury. Sci Rep. 2016;6:29121. https://doi.org/10.1038/srep29121.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Li S, Wang X, Gu Y, Chen C, Wang Y, Liu J, Hu W, Yu B, Wang Y, Ding F, Liu Y, Gu X. Let-7 microRNAs regenerate peripheral nerve regeneration by targeting nerve growth factor. Mol Ther. 2015;23(3):423–33. https://doi.org/10.1038/mt.2014.220.

    Article  CAS  PubMed  Google Scholar 

  34. Yi S, Wang QH, Zhao LL, Qin J, Wang YX, Yu B, Zhou SL. miR-30c promotes Schwann cell remyelination following peripheral nerve injury. Neural Regen Res. 2017;12(10):1708–15. https://doi.org/10.4103/1673-5374.217351.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011;146(3):353–8. https://doi.org/10.1016/j.cell.2011.07.014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Tay Y, Kats L, Salmena L, Weiss D, Tan SM, Ala U, Karreth F, Poliseno L, Provero P, Di Cunto F, Lieberman J, Rigoutsos I, Pandolfi PP. Coding-independent regulation of the tumor suppressor PTEN by competing endogenous mRNAs. Cell. 2011;147(2):344–57. https://doi.org/10.1016/j.cell.2011.09.029.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Karreth FA, Tay Y, Perna D, Ala U, Tan SM, Rust AG, DeNicola G, Webster KA, Weiss D, Perez-Mancera PA, Krauthammer M, Halaban R, Provero P, Adams DJ, Tuveson DA, Pandolfi PP. In vivo identification of tumor- suppressive PTEN ceRNAs in an oncogenic BRAF-induced mouse model of melanoma. Cell. 2011;147(2):382–95. https://doi.org/10.1016/j.cell.2011.09.032.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, Tramontano A, Bozzoni I. A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell. 2011;147(2):358–69. https://doi.org/10.1016/j.cell.2011.09.028.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Morey JS, Ryan JC, Van Dolah FM. Microarray validation: factors influencing correlation between oligonucleotide microarrays and real-time PCR. Biol Proced Online. 2006;8:175–93. https://doi.org/10.1251/bpo126.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Beckman KB, Lee KY, Golden T, Melov S. Gene expression profiling in mitochondrial disease: assessment of microarray accuracy by high-throughput Q-PCR. Mitochondrion. 2004;4(5–6):453–70. https://doi.org/10.1016/j.mito.2004.07.029.

    Article  CAS  PubMed  Google Scholar 

  41. Etienne W, Meyer MH, Peppers J, Meyer RA Jr. Comparison of mRNA gene expression by RT-PCR and DNA microarray. Biotechniques. 2004;36(4):618–20, 622, 624–616. https://doi.org/10.2144/04364ST02.

    Article  CAS  PubMed  Google Scholar 

  42. White TE, Surles-Zeigler MC, Ford GD, Gates AS, Davids B, Distel T, LaPlaca MC, Ford BD. Bilateral gene interaction hierarchy analysis of the cell death gene response emphasizes the significance of cell cycle genes following unilateral traumatic brain injury. BMC Genomics. 2016;17:130. https://doi.org/10.1186/s12864-016-2412-0.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by the National Natural Science Foundation of China [31700926], and the Priority Academic Program Development of Jiangsu Higher Education Institutions of China [PAPD].

Availability of data and materials

Sequencing data were uploaded to NCBI database (accession number PRJNA394957; SRP113121).

Author information

Author notes

  1. Tianmei Qian and Chunlin Fan contributed equally to this work.

Authors and Affiliations

  1. Key Laboratory of Neuroregeneration of Jiangsu and Ministry of Education, Co-innovation Center of Neuroregeneration, Nantong University, 19 Qixiu Road, Nantong, Jiangsu Province, People’s Republic of China, 226001

    Tianmei Qian, Qianyan Liu & Sheng Yi

  2. Department of Respiratory and Critical Care Medicine, Affiliated Hospital of Nantong University, Nantong, Jiangsu, China

    Chunlin Fan

Authors
  1. Tianmei Qian
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Chunlin Fan
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Qianyan Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Sheng Yi
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceived and designed the experiments: SY. Performed the experiments: TQ QL SY. Analyzed the data: CF SY. Contributed reagents/materials/analysis tools: SY. Wrote the manuscript: CF SY. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Sheng Yi.

Ethics declarations

Ethics approval

All the experimental procedures involving animals were conducted in accordance with Institutional Animal Care guidelines of Nantong University and approved ethically by the Administration Committee of Experimental Animals, Jiangsu, China.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. List of primer pairs for RT-PCR. (XLSX 19 kb)

Additional file 2:

Table S2. List of significantly differentially expressed RNAs in the sciatic nerve segments after injury. RNAs with fold change> 10 or < − 10 and FDR < 0.001 as compared with 0 day control were considered as significantly differentially expressed. Gene ID, gene symbol, log2Ratio, fold change, p-value, and FDR of significantly differentially expressed mRNAs and lncRNAs were listed. (XLSX 379 kb)

Additional file 3:

Table S3. List of enriched GO terms of significantly differentially expressed mRNAs in the sciatic nerve segments after injury. (XLSX 124 kb)

Additional file 4:

Table S4. List of enriched GO terms and KEGG pathways of co-expressed mRNAs and lncRNAs. (XLSX 86 kb)

Additional file 5:

Table S5. List of enriched GO terms of co-expressed mRNAs and lncRNAs with validated mRNAs. (XLSX 23 kb)

Additional file 6:

Table S6. RPKM reads of mRNAs and lncRNAs in the sciatic nerve segments after injury. Expression levels of mRNAs (HMOX1 and LIF) and lncRNAs (XLOC_138166, XLOC_150234, XLOC_083376, XLOC_174287, XLOC_146853, XLOC_065278, XLOC_130885, XLOC_000412, XLOC_083385, XLOC_027762, XLOC_146842, XLOC_174535, XLOC_108671, XLOC_174536, XLOC_172336, XLOC_056487, XLOC_134667, XLOC_133837, XLOC_150864, and XLOC_080446) in the sciatic nerve segments at 0, 1, 4, 7, and 14 days after injury were demonstrated as RPKM reads. (XLSX 12 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Qian, T., Fan, C., Liu, Q. et al. Systemic functional enrichment and ceRNA network identification following peripheral nerve injury. Mol Brain 11, 73 (2018). https://doi.org/10.1186/s13041-018-0421-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13041-018-0421-4

Keywords

Molecular Brain

ISSN: 1756-6606