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

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 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 cis-target 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.

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 post-nerve 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 (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 unique DElncRNAs and common DElncRNAs were both enriched in the PI3K-Akt signaling pathway. However, the specific DElncRNAs were totally different in these two groups: the control vs D4 unique DElncRNAs that were enriched in the PI3K-Akt signaling pathway included Col4a5, Fgf1, Ghr, Magi2, Met, Pdgfa, Ppp2r2b, Ppp2r2c, Prkaa2, Prlr and Thbs2, while the common DElncRNAs enriched in this pathway included Ccne1, Cdkn1a, Col9a1, Erbb4, Fgf7, Itgb8, Lamc2 and 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 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) most lncRNAs are consistent with their corresponding mRNA regulation trends, and the expression trends of only a few pairs of lncRNAs and mRNAs are opposite, including NONRATT028679.2-Smad3, NONRATT028680.2-Smad3, NONRATT015315.2-Igsf10 and NONRATT030680.2-Col4a5; (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.