Integrated analysis of long noncoding RNAs and mRNA expression profiles reveals the potential role of lncRNAs in early stage of post-peripheral nerve injury in Sprague-Dawley rats

The regulatory role of lncRNAs in the early stage post peripheral nerve injury (PNI) is not yet clear. In this study, next-generation sequencing was used to perform deep sequencing on normal sciatic nerves (control) and lesional tissues derived on the 4th (D4) and 7th days (D7) after sciatic nerve injury in rats. Time-point unique differentially expressed lncRNAs (DElncRNAs) were analyzed for functional enrichment. The results showed that 776 DElncRNAs were unique to D4, and their functions were mainly enriched in wound healing, phosphatase binding and MAPK signaling pathways; 317 DElncRNAs were unique to D7, and their functions were mainly enriched in ion transmembrane transporter channel activity; 579 DElncRNAs were shared by these two days, and their functions were mainly enriched in axongenesis, the PI3K-Akt signaling pathway and cell cycle. Furthermore, lncRNA-mRNA interaction networks were constructed in functions or pathways with a high enrichment rate. Finally, 3 mRNAs and 4 lncRNAs in the axongenesis interaction network were selected, and their expression levels were verified by RT-qPCR. This study preliminarily revealed the regulatory role of lncRNAs at different time points in the early stage post PNI, which provides potential targets for basic research and clinical treatment of PNI.


INTRODUCTION
The morbidity of peripheral nerve injury (PNI) has been increasing with the development of social industrialization, and such an injury has a negative effect on the individual and family life, causes the loss of social labor, and requires a large amount of medical resources [1,2]. The main clinical manifestations of the disease include peripheral sensory, motor and autonomic nerve disorders [3]. The pathophysiological changes after PNI include Wallerian degeneration of the distal stump of the injured nerve and the formation of a growth cone from the proximal stump, followed by the regenerated proximal nerve fibers growing through to the distal end and reinnervating the target organs to complete the nerve regeneration process [4]. The expression of a number of pivotal genes, enriched in the inflammatory response, immune response, cell migration, cell proliferation and other key biological processes involved in peripheral nerve regeneration was reported to reach a high plateau within one week [5]. During this regeneration and repair process, such genes undergo a complex and dynamic regulation [6]; however, to date, few studies have examined this issue. Therefore, a comprehensive understanding of the regulation rules in such a process could provide a AGING solid foundation for the basic research and clinical application of PNI treatment.
LncRNAs are a type of noncoding RNA with a length of more than 200 nucleotides [7]. Under usual circumstances, lncRNAs do not encode proteins themselves, but they can regulate gene expression in different ways, and in turn participate in many physiological and pathological processes [8]. According to the mechanism of gene expression regulation, lncRNAs can be roughly divided into two categories: cis and trans. Cis-acting lncRNAs mainly regulate their neighboring chromatin state or gene expression, while trans-acting lncRNAs affect gene expression in remote transcripts [9].
Limited research on the regulatory role of lncRNAs in the process of peripheral nerve injury and repair has been reported thus far. Ma and colleagues found that the expression of lncRNA MEG3 was upregulated in a sciatic nerve injury model, and downregulating the expression of such lncRNA in Schwann cells can promote cell proliferation and migration and induce the growth of nerve axons [10]. Another study showed that lncRNA TNXA-PS1 could function as an endogenous competing RNA, which regulates the expression of Dusp1 protein in Schwann cells by adsorbing miR-24-3p and miR-152-3p and in turn affects the migration of Schwann cells [11]. Pan and other researchers used microarray data to analyze the regulatory relationship between differentially expressed lncRNAs and mRNAs in the early stage of mouse sciatic nerve injury. The results of the analysis suggest that JUN, a key regulatory gene after PNI, could be cis-regulated by lncRNA ENSMUSG00000087366, and the other gene MBP, mainly involved in myelination, may be transregulated by lncRNA ENSMUSG00000084785 [12]. However, most of the abovementioned studies were based on gene chip technology. The main disadvantage of such a technology is that it can only detect and analyze known genes and omits unknown lncRNAs that may play a role in peripheral nerve damage and repair.
In this study, we used next-generation sequencing (NGS) to perform deep sequencing of lncRNAs and mRNAs at different time points after rat sciatic nerve injury. Differentially expressed lncRNAs (DElncRNAs), unique to each time point or shared by both, were first analyzed and combined with the corresponding differentially expressed mRNA (DEG) data we predicted the cis-and trans-acting target genes of the DElncRNAs. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were then performed on these DElncRNAs, and we constructed lncRNA-mRNA interaction networks in some GO terms and KEGG pathways with a high enrichment score. Finally, the expression of 3 mRNAs and 4 lncRNAs in network axongenesis was verified by RT-qPCR. This study preliminarily revealed the potential regulatory role of lncRNAs at different time points in the early stage of PNI.

Laboratory animal selection and ethics statement
Twenty-one 1-month-old Sprague-Dawley (SD) rats weighing approximately 150-200 g were randomly divided into 3 groups: the control group (control), 4 days after sciatic nerve injury (D4) and 7 days after injury (D7). This experiment was approved by the Experimental Animal Ethics Committee of Sun Yat-sen University. All experimental operations complied with the rules and regulations formulated by the ethics committee, and the damage and pain caused to animals were minimized during the experiment. The overall experimental design and workflow of this research are shown in Figure 1.

Sciatic nerve clamp injury model
Animals in the experimental groups (D4 and D7) were anesthetized by intraperitoneal injection of 10% chloral hydrate (0.3 ml/100 g) and were subjected to sciatic nerve crush injury model according to our previously published articles [13,14]. The operation is briefly described as follows. After the left hind limb and back of the rat were shaved and disinfected, the skin was cut parallel to the femur, and the superficial gluteus muscle was bluntly separated to fully expose the left sciatic nerve. The middle segment of the sciatic nerve was clamped with hemostatic forceps for 30 s and then released immediately, and the wound was closed layer by layer. Animals in the control group underwent the same operation but no sciatic nerve clamp operation. All rats were kept in a suitable and stable environment, ambient temperature between 20-25° C and humidity between 50-65%, with adequate food and water supply.

RNA extraction and purification
Two rats were randomly selected from each group, and RNA extraction and purification were performed on either normal or injured sciatic nerve samples using the mir VanaTM miRNA Isolation Kit (Cat# AM1561, Ambion, MA, USA) according to the standard protocol. The total RNA obtained was subjected to quality inspection by electrophoresis and purification using an RNAClean XP kit (Cat# A63987, Beckman Coulter, Inc., CA, USA) and RNase-Free DNase Set (Cat#79254, QIAGEN, Duesseldorf, Germany). The purified RNA was tested for quality by NanoDrop ND-2000 spectrophotometer (Thermo Fisher, MA, USA) and AGING Agilent Bioanalyzer 2100 (Agilent Technologies, CA, USA).

Sequencing library construction, RNA sequencing and data analysis
Purified total RNA was subjected to rRNA removal, fragmentation, first-and second-strand cDNA synthesis, end repair, poly-A tail addition, ligation, enrichment and other steps to construct a sequencing sample library. High-throughput RNA paired-end sequencing was then performed using the Illumina HiSeq Xten platform (Illumina, Inc., CA, USA). The unqualified reads, with low overall or end quality and containing sequencing primers, etc. were removed from the raw reads to obtain clean reads, from which rat ribosomal RNAs were further filtered. The spliced mapping algorithm in Hisat2 (version: 2.0.4) [15] was applied to map the above processed reads to the rat reference genome (ftp://ftp.ensembl.org/pub/release-100/fasta/ rattus_norvegicus/dna/).
With the use of the Stringtie software (version: 1.3.0) [16], the mapped reads were assembled to construct the transcripts, which were then matched to the reference databases (NOCODE and Ensembl) using GffCompare (version: 0.9.8) software to obtain annotation. The unmatched transcripts were subjected to potential novel lncRNA prediction with the following criteria: (1) transcript length ≥ 200 bp AND exon ≥ 2, (2) the predicted ORF < 300 bp, (3) coding potential calculator (CPC) score < 0 [17] AND coding-noncoding index (CNCI) score < 0 [18] AND with no significant difference using Pfam comparison [19]. The Perl script was then applied to locate the position of the above novel lncRNAs on the chromosomes and to make annotations.
With the use of Stringtie software, the fragment numbers of each gene were counted, which were then normalized by TMM (trimmed mean of M values) method. Finally the FPKM value of each known or newly predicted lncRNAs was calculated using perl script, where the lncRNA ID starting with NON is the known lncRNA in the NONCODE database, the ID starting with ENS is the known lncRNA in the Ensembl database, and the ID starting with MSTRG is the newly predicted lncRNA. The edgeR package [20] was then used to analyze the differentially expressed genes between samples (Control vs D4 and Control vs D7). The screening criteria for differentially expressed genes were as follows: (1) q-value (the corrected p-value) < 0.05, and (2) the absolute value of fold change ≥ 2. The RNA-seq data were released to the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database under the accession number GSE162548.

Differentially expressed lncRNA grouping and target gene prediction
The intersection of the DElncRNAs, compared between the control vs D4 and the control vs D7 groups, was made, and the DElncRNAs were divided into three groups according to the intersection results: (1) DElncRNAs unique to the control vs D4 group, (2) common DElncRNAs shared by the control vs D4 and the control vs D7 groups, and (3) DElncRNAs unique to the control vs D7 group. Cis-and trans-target gene prediction was then performed for the DElncRNAs in the above three groups. Target genes located with a distance less than 10 kb from the upstream or downstream of the lncRNA were considered as the cistarget gene. For trans-target gene prediction, the NCBI BLAST tool was first used to select rat mRNA sequences that had complementarity or similarity with the DElncRNA sequences, and then the RNAplex software [21] was used to calculate the complementary energy between the two sequences. Target sequences with an identity greater than 85% and an e-value greater than e-20 were selected as the trans-targets of the lncRNA. Refer to the Supplementary Tables 1, 2 for the cis-and trans-target genes, respectively.

LncRNA functional enrichment and lncRNA-mRNA interaction network construction
According to the DElncRNAs grouping method, the differentially expressed mRNAs (DEGs, DEG screening criteria were the same as DElncRNAs), compared between the control vs D4 and the control vs D7 groups, were intersected and then divided into three groups: (1) DEGs unique to the control vs D4 group, (2) common DEGs shared by these two groups, (3) DEGs unique to the control vs D7 group. The predicted target genes of DElncRNAs intersecting with the corresponding DEGs in each group were subjected to GO functional enrichment analysis, including biological process (BP), molecular function (MF) and cellular component (CC), and KEGG pathway enrichment analysis using the Bioconductor program [22] in R software. The enrichment results with p-values less than 0.05 were selected and displayed with dot plots. The GO terms or KEGG pathways with a high enrichment score in each group were selected for lncRNA-mRNA interaction network analysis, and the results were processed and displayed with Cytoscape software.

Real-time quantitative PCR verification
Three sets of lncRNA-mRNA interaction networks in the GO term axongenesis were selected, namely, Gndf and NONRATT015075.2, Nfasc and NONRATT008698.2, Pmp22 and NONRATT004387.2, NONRATT004386.2. The gene expression of the above 3 mRNAs and 4 lncRNAs was verified by RT-qPCR in the remaining 5 animals in each group. Reverse transcription of RNA was performed using the ReverTra Ace qPCR Kit (TOYOBO, #FSQ-101, Osaka, Japan). The reaction conditions were set as follows: 37° C, 15 min, then 98° C, 5 min, and the obtained cDNA was stored at 4° C. The PCR system was prepared according to the Power SYBR Green PCR Master Mix (ABI, #4368708, MA, USA) kit operation manual and then subjected to RT-qPCR in the QuantStudio 5 Real-Time PCR System (ABI) instrument. The reaction conditions were set as follows: 50° C, 2 min, then 95° C, 10 min, and then 40 reaction cycles of 95° C, 15 s and 60° C, 1 min. Refer to Supplementary Table 3 for the PCR primer sequences.

Statistical analysis
The RT-qPCR data were statistically analyzed by R software, and the statistical results are expressed as the mean ± standard error (SEM). If the data of two independent samples were content with normal distribution and the variances were uniform, then Student's t-test was used to analyze the differences between these two samples. If the variances were not uniform, the t-test with Welch's correction was used for analysis. If the data were not normally distributed, the Wilcoxon test was used for analysis. The results were displayed using R software and GraphPad Prism software (version: 8.4.3). A p-value less than 0.05 was considered statistically significant, such as *P < 0.05, **P < 0.01, ***P < 0.001. The RT-qPCR statistical results were all from 5 independent experiments.

RESULTS
After completing the construction of the cDNA library of the sequencing samples, RNA sequencing was performed using the Illumina HiSeq Xten platform. The number of raw reads in each group and subsequent processed data are shown in Table 1. Some unqualified reads with low overall quality, sequencing primers, low end quality, etc., were removed to obtain clean reads. The results showed that the clean reads accounted for a high proportion of the raw reads (94.76%-96.43%), indicating the high quality and reliability of our sequencing data. Ribosomal RNA reads and reads with only single-end sequencing results were further removed, and the remaining reads (pre-mapped reads) were matched to the rat genome sequences to obtain the mapped reads. The mapping rate of all the 6 groups reached a percentage greater than 90%, except for the Con_a group (88.70%). The unique mapped reads ratio of each group was above 99%.

LncRNA characteristics analysis
The lncRNAs identified by sequencing were classified and counted according to their lengths ( Figure 2A). As the length of lncRNAs increases, its proportion to the total lncRNA number decreases. LncRNAs with a length of 200 -1,000 nt accounted for the highest proportion, up to 75.03%, and lncRNAs with a length greater than 10,000 nt accounted for the lowest proportion, only 0.97%. Furthermore, lncRNAs were classified according to their location in the genome ( Figure 2B), and the results showed that the lncRNAs originating from the exon region of the coding sequences (exonic sense) accounted for the highest proportion (54.6%), followed by the intergenic lncRNAs (21.8%). The lncRNAs originating from the intron region of the antisense strand (intronic antisense) accounted for the lowest proportion (1.8%). The lncRNAs and mRNAs were then statistically analyzed according to the exon numbers ( Figure 2C), and the results showed that lncRNAs containing only one exon accounted for the largest proportion (38.86%), and the proportion of lncRNAs containing 1 to 4 exon numbers to the total number of lncRNAs was over 91%. Such results were consistent with the lncRNA length proportion: the proportion of lncRNAs less than 2,000 nt was 90.41%. Finally, the average expression value of lncRNAs and mRNAs was analyzed ( Figure 2D): The average expression level of lncRNAs is slightly higher than that of mRNAs, and the expression value distribution is more concentrated.

DElncRNAs analysis and grouping
The expression values of lncRNAs on D4 and D7 postnerve injury were compared to the normal tissue (control vs D4 and control vs D7). In the control vs D4 group, a total of 1355 DElncRNAs were identified, of which 691 were upregulated and 664 were downregulated ( Figure 3A), and in the control vs D7 group, a total of 896 DElncRNAs were identified, of which 431 were upregulated and 465 were downregulated ( Figure 3B). The same screening criteria were applied to identify the DEGs: in the control vs D4 group, a total of 3607 DEGs were identified, of which 2240 were upregulated and 1367 were downregulated ( Figure 3C), and in the control vs D7 group, a total of 2990 DEGs were identified, of which 1798 were upregulated and 1192 were downregulated ( Figure 3D).
To further investigate the role that lncRNAs play at different time points in the early stage of post-sciatic nerve injury, we intersected the DElncRNAs of control vs D4 and control vs D7 groups and identified some differences and similarities: there are 776 DElncRNAs unique to the control vs D4 group, 317 DElncRNAs unique to the control vs D7 group, and 579 common DElncRNAs shared by these two groups. The obtained data were sorted to prepare for the subsequent functional enrichment analysis ( Figure 3E), and the unique and common DEGs were analyzed and sorted using the same method ( Figure 3F).

GO and KEGG pathway enrichment of the DElncRNAs
The cis-and trans-target genes (Supplementary Materials) were identified on the DElncRNAs unique to the control vs D4 group, the DElncRNAs unique to the control vs D7 group, and the common DElncRNAs in these two groups. After intersecting the DEGs in the corresponding groups, these target genes were subjected to GO and KEGG pathway enrichment analysis to better understand the function of the DElncRNAs in each group. The screening criteria were set as follows: the enrichment p value was less than 0.05 (note that the legends in Figures 4, 5 display the adjust p value).
The results of GO enrichment analysis showed that in BP the DElncRNAs unique to the control vs D4 were mainly enriched in wound healing, regulation of response to wound and endothelial cell chemotaxis and other GO terms related to wound healing; in CC, these lncRNAs were mainly enriched in cellular regions related to chromatin, centriole and sarcoplasmic reticulum; and they were only enriched in a single GO term, phosphatase binding, in MF ( Figure 4A and Table 2). The common DElncRNAs shared by control vs D4 and control vs D7 groups were mainly enriched in axongenesis, chromosome segregation and regeneration and other biological processes related to axon regeneration and cell cycle; in CC, they were mainly enriched in spindle, presynaptic membrane and main axon and other regions; in MF, they were enriched in actin binding, heparin binding and sodium and potassium ion channel activity AGING Clean reads ratio = Clean reads/Raw reads; Mapped reads ratio = Mapped reads/Pre-mapped reads; Pre-mapped reads = Clean reads -rRNA reads -single-end reads; Unique mapped reads ratio = unique mapped reads/mapped reads.
( Figure 4B and Table 3). The DElncRNAs unique to the control vs D7 group were mainly enriched in ion transmembrane transporter channel activity in BP; these lncRNAs were mainly enriched in transmembrane transport complex and ion channel complex in CC; and they were mainly enriched in ion channel binding and potassium ion transmembrane transporter activity in MF ( Figure 4C and Table 4).
KEGG pathway enrichment results (Table 5) showed that the DElncRNAs unique to the control vs D4 group were mainly enriched in the PI3K-Akt signaling pathway, MAPK signaling pathway and Rap1 signaling pathway ( Figure 5A). The common DElncRNAs of these two groups were mainly enriched in the PI3K-Akt signaling pathway, focal adhesion and the cell cycle ( Figure 5B). Interestingly, the control vs D4 group  Tnr. The DElncRNAs unique to the control vs D7 group were mainly enriched in the long-term potentiation, hematopoietic cell lineage and aldosterone synthesis and secretion pathways ( Figure 5C).

lncRNA-mRNA interaction network construction
The GO terms and KEGG pathways with high enrichment scores and strong correlations with nerve  AGING regeneration were selected, and lncRNA-mRNA interaction networks were constructed within each group, with the regulation trend and the interaction mechanism annotated. The lncRNA-mRNA interaction networks of wound healing and the PI3K-Akt pathway, highly enriched GO term and KEGG pathway of the control vs D4 group unique DElncRNAs are displayed in Figure 6A, 6B, respectively. The networks of axongenesis and cell cycle pathways, highly enriched GO term and KEGG pathway of common DElncRNAs are displayed in Figure 6C, 6D, respectively. The networks of ion transportation and hematopoiesis, highly enriched GO term and KEGG pathway of the control vs D7 group unique DElncRNAs are displayed in Figure 6E, 6F, respectively. From the interaction networks, we found the following rules: (1)  (2) the action mode of most lncRNAs on mRNAs is cis-acting, except NONRATT015315.2-Igsf10 and NONRATT026281.2-Rapgef3, whose mode of action is trans-acting; (3) the regulatory relationship between most lncRNAs and mRNAs is one-to-one, only a few mRNAs are regulated by two lncRNAs, and very few mRNAs are regulated by three lncRNAs at the same time.

RT-qPCR verification
Three sets of lncRNA-mRNA interaction pairs enriched in GO term axongenesis ( Figure 6C), including 3 mRNAs and 4 lncRNAs, namely, Gndf and NONRATT015075.2, Nfasc and NONRATT008698.2, Pmp22 and NONRATT004387.2, NONRATT004386.2 were selected, and their expression was measured using RT-qPCR to verify the accuracy of the RNA-seq sequencing results. The RT-qPCR results (Figure 7) showed that the expression trends of the above 7 genes on D4 and D7 post-operation compared to the control group were completely consistent with the RNA-seq results. In addition to the NONRATT008698.2 gene, the expression levels of the other 6 genes on D4 and D7 all had significant differences compared to the control group.

DISCUSSION
After peripheral nerve injury, the injured distal axon undergoes Wallerian degeneration, accompanied by a phenotypic change in Schwann cells from dedifferentiation (to guide axon regeneration) to redifferentiation (to promote the myelination of regenerated axons) [23], and this process involves a variety of biological activities from gene transcription and epigenetic regulation [24,25]. Although an increasing number of studies have attempted to clarify the role of a single gene in this process in recent years, the regulatory relationship among genes is still not clear. In this study, with the use of RNA-seq techniques and bioinformatics methods, we identified the differentially expressed lncRNAs, predicted their potential role in the early post-PNI stage, and finally constructed the lncRNA-mRNA interaction networks in the pathways or biological processes in which those lncRNAs were enriched.
In our study, 579 common DElncRNAs were shared by the control vs D4 and control vs D7 groups, which participate in the whole process of the regulation of gene expression in the early stage after sciatic nerve injury. These lncRNAs were mainly enriched in axongenesis, chromosome segregation, regulation of cell development, spindle, presynaptic membrane and ion channel activity in GO enrichment. In KEGG enrichment analysis, these DElncRNAs were mainly enriched in the PI3K-Akt signaling pathway, focal adhesion, JAK-STAT signaling pathway, the ECM-receptor interaction pathways. These     GO terms or KEGG pathways that DElncRNAs were mainly enriched were highly consistent with the results of the DEG functional enrichment analysis in rat sciatic nerve injury reported by previous studies [26][27][28], indicating a high credibility of evidence that the lncRNAs we analyzed have the potential to play a role in the regulation of PNI-associated genes. To determine the different functions of the DElncRNAs at different time points after injury, we analyzed the GO and KEGG enrichment of DElncRNAs unique to the control vs D4 and control vs D7 groups. The results showed that the unique DElncRNAs in the control vs D4 group were mainly enriched in wound healing, regulation of wound healing and endothelial cell chemotaxis-related GO terms, suggesting that the regulatory effect of lncRNAs on injured nerve repair and angiogenesis may begin in the very early stages after nerve injury. The DElncRNAs unique to the control vs D7 group were mainly enriched in GO terms related to ion transmembrane transport, especially channels related to calcium and potassium, which play important roles in axon degeneration and nerve repair. After PNI, the influx of extracellular calcium and the release of intracellular calcium can induce the activation of axon proteases, such as calpains, which in turn cause Wallerian degeneration at the distal end of the injury [29]. Studies have shown that the reduction in intracellular calcium can lead to the delayed occurrence of axon degeneration [30]. Studies on AGING potassium channels in PNI have shown that Kv2.1 ion channels on the surface of motor neurons rapidly decrease one week after injury and then gradually recover [31], although the influence of such a change is not clear. In the GO term axonogenesis enriched by the common DElncRNAs, lncRNA NONRATT015075.2 may lead to the upregulation of the glial cell-line-derived growth factor (Gdnf) gene. Gdnf, a TGF-β-related neurotrophic factor, can promote the proliferation and migration of Schwann cells and help the regeneration of neuronal axons [36][37][38]. Sema3e (semaphorin 3e) and its receptor play an important role in axonal guidance [39]. In this study, we predicted 3 lncRNAs that could regulate Sema3e, namely, NONRATT021162.  [40]. Finally, we selected 3 groups of lncRNA-mRNA interaction pairs, a total of 7 genes including 3 mRNAs and 4 lncRNAs, and used RT-qPCR to detect their expression levels. The results showed that the results of RT-qPCR and RNA-seq results were completely consistent. However, the regulatory relationship between lncRNAs and their target genes needs further study.
Studies have shown that aging and its comorbidities, such as hypertension and diabetes, have an important impact on the plasticity of both central and peripheral nerve systems after injury [41][42][43]. However, in the current study, as we did not include the aging animals as our study objects, it is unable for us to clarify the impact of senescence on the lncRNAs expression profile after peripheral nerve injury, which is the limitation of this study and needs us to further investigate in the future studies.
In this study, we used NGS technology to sequence rat sciatic nerve tissues on D4 and D7 after injury, analyzed the DElncRNAs at different time points, and predicted the target genes combined with the corresponding mRNA expression data. The DElncRNAs at different time points were analyzed for functional enrichment, and the GO terms and KEGG pathways with high enrichment rates were selected to construct the lncRNA-mRNA interaction networks. This study preliminarily revealed the regulatory role of lncRNAs at different time points in the early stage after PNI, which could provide new potential targets for the research and treatment of PNI. were selected from the lncRNA-mRNA network constructed in axongenesis GO term, and their expression, expressed as fold change value compared with control group (omitted), were verified using RT-qPCR (red bar). The lncRNA-mRNA pairs selected are Gdnf (A) and NONRATT015075.2 (D); Nfasc (B) and NONRATT008698.2 (E); Pmp22 (C), NONRATT004386.2 (F) and NONRATT004387.2 (G). And the results were compared with those obtained from RNA-seq (blue bar). Statistically significant values are presented as asterisks (*), *P < 0.05, **P < 0.01.