- Open Access
Systemic functional enrichment and ceRNA network identification following peripheral nerve injury
Molecular Brain volume 11, Article number: 73 (2018)
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.
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 . 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 . 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 . 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 .
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 . 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.
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.
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 . 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.
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 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.
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 . 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).
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.
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.
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).
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).
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).
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.
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 . 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.
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 . 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 . 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 . 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.
One-way analysis of variance
Competing endogenous RNA
Database for Annotation, Visualization, and Integrated Discovery
False discover rate
Heme oxygenase (decycling) 1
Kyoto Enrichment of Genes and Genomes
Leukemia inhibitory factor
Long non-coding RNA
Peripheral myelin protein 22
Reads per kilobase transcriptome per million mapped reads
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Oliveros JC. VENNY. An interactive tool for comparing lists with Venn Diagrams; 2007.
Guttman M, Rinn JL. Modular regulatory principles of large non-coding RNAs. Nature. 2012;482(7385):339–46. https://doi.org/10.1038/nature10887.
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.
Kurokawa R, Rosenfeld MG, Glass CK. Transcriptional regulation through noncoding RNAs and epigenetic modifications. RNA Biol. 2009;6(3):233–6.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
This work was supported by the National Natural Science Foundation of China , 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).
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
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Table S1. List of primer pairs for RT-PCR. (XLSX 19 kb)
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)
Table S3. List of enriched GO terms of significantly differentially expressed mRNAs in the sciatic nerve segments after injury. (XLSX 124 kb)
Table S4. List of enriched GO terms and KEGG pathways of co-expressed mRNAs and lncRNAs. (XLSX 86 kb)
Table S5. List of enriched GO terms of co-expressed mRNAs and lncRNAs with validated mRNAs. (XLSX 23 kb)
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)
About this article
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
- Peripheral nerve injury
- RNA sequencing
- Bioinformatic analysis