The raw microarray data were downloaded from the Gene Expression Omnibus database¹ (GEO), platform Agilent-014879: Whole Rat Genome Microarray 4 × 44 K G4131F. Gene expression signatures for 30 samples were generated, including L4-6 DRG and proximal SN tissues (0.5 cm) at day 0, day 1, day 4, day 7, and day 14 after SN resection. This dataset included three samples in each tissue, DRG, and SN, and the expression levels were measured for all samples at five time points.

The differential expression was computed with GeneSpring (GX V 12.6). The time baseline was assigned the same way in both studies, with DRG and SN. Given the start day, 0 day = b, further measurements that were taken at successive intervals were redefined as: I1 = 1 day − b, I2 = 4 day − b, I3 = 7 day − b, I4 = 14 day − b. Sample groups thus appear according to time point in each tissue, for a total of five groups. Significance was assigned when differences in expression values between the groups had fold change (FC) >1.5, and t-test p-values <0.05, then adjusted for multiple testing by the false discovery rate (FDR) or Benjamini and Hochberg method (24). The probe sets and genes correspondence was established by FDR cutoff q < 0.05, and this is a standard way to bypass the problem of probes mapping to multiple genes. The outcome of this process (see also Table S8 in Supplementary Material) was the gene set called DEGs, retained for downstream analysis.

Of course, the stringency of these criteria influences the relative abundance or paucity of evidences that allow downstream analysis and thus inference; however, the goal of capturing biotypes other than genes, say ncRNA, suggests considering that relatively low expression values should be expected to characterize them.

Pathway analysis was obtained using the tool GeneGo Metacore™ (Thomson Reuters Corporation, New York, NY, USA). The DEGs identified by GeneSpring analysis were uploaded to the GeneGo website, and core analysis was applied to obtain the list of activated pathways, using 0.05 as a threshold for p-values. While we rely on established tools for protein-coding gene annotations, i.e., DEGs enriching for pathway terms, we need to work at a more empirical level for the ncRNAs, as explained below.

The differentially expressed detections were annotated with both protein-coding and non-coding biotypes, as provided by Ensembl rel. 77 (25), Rattus_norvegicus.Rnor_6.0.79.gtf, and Mus_musculus.GRCm38.76.gtf. In the non-coding category, annotations regarding lncRNAs were provided for pseudogenes, lincRNAs, and asRNA – considered of major interest in our work. Due to the categories appearing under the label pseudogene, we aimed to have one single group, and thus the processed pseudogenes were consolidated into the final pseudogene group. Below, we explain the annotation strategy for these three major categories. Lastly, we indicate the use of targets to build protein–protein interaction networks (PINs).

The parental protein genes found differentially expressed in either of the neuronal tissues were used to generate the PIN using STRING² (V. 9.1). This is a known state-of-the-art database reporting both experimental and predicted protein–protein interactions. This network was functional to our downstream contextual analysis, i.e., the contextualization of the identified targets with respect to biological processes and pathways.

We used Mus_musculus.GRCm38.cdna.all.fa and Rattus_norvegicus.Rnor_6.0.cdna.all.fa from Ensembl (rel 77). The mouse was used as a reference only.

The microarray dataset was profiled for the DEGs in both neuronal tissues – DRG and SN – at different time points (see Table 1). In addition, annotated differentially expressed probes and ncRNAs are uniquely obtained using the mouse reference (see Materials and Methods for details and bold numbers in parentheses reported in Table 1).

We detected 32 unique lincRNAs in SN and 8 unique lincRNAs in DRG. Then, we detected 31 unique antisense in SN and 12 unique antisense in DRG. Such antisense detections in SN and in DRG were not significantly associated with targets; hence, they were no longer explored. We additionally identified 452 unique pseudogenes in SN and 56 unique pseudogenes in DRG, subsequently investigated for their protein-coding targets.

Pathway analysis was conducted with the tool GeneGo, delivering evidences displayed in Figures 1A,B and 2A,B and summarized in the top-10 annotated terms for both tissues, DRG and SN. The enriching DEG sets appear in supplemental files (see Tables S6 and S7 in Supplementary Material for complete lists of annotated terms). Also, fractions of genes enriching pathway terms are reported inside the plots, while the listed terms are sorted by FDR-corrected enrichment values (from best, at the top).

We can observe among the terms in Figure 1A the one annotated as neurophysiological process_Kappa-type opioid receptor in transmission of nerve impulses, especially involved in DRG_I2. The relevance of this term is associated to CREM (activators), AP-1, and P/Q-type calcium channel alpha-1A subunit. The latter gene is a voltage-dependent calcium channel that mediates the entry of calcium ions into excitable cells and is involved in various calcium-dependent processes, like muscle contraction, neurotransmitter release, and gene expression. It is also primarily expressed in neuronal tissue, and it was shown that mutations to this gene could potentially cause two neurological disorders (29).

We also noted recurrence of cell adhesion_Integrin inside-out signaling in T cells at SN_I2, SN_I3, and SN_4, but not at SN_I1. This may indicate a possible role for the integrin-regulated cell adhesion in tissue morphogenesis and wound healing, together with the regulation of cell growth and differentiation, whereas the temporal pattern points to recurrence after injury.

Another interesting pathway is immune_response_CCL2_signaling, which is involved mainly in SN and refers to chemokines previously shown to be upregulated in injured DRG neurons under different NP models (30–33). We found enrichment of immune response_NFAT in SN_I2, SN_I4 and signal transduction_cAMP signaling in DRG_I4 by calcineurin A (catalytic), which is a calcium and calmodulin dependent serine/threonine protein phosphatase. This gene is relevant for CNS functions, such as neurite extension, synaptic plasticity, learning, and memory (34, 35).

Most pathways in SN show gene enrichment for ERK1/2, a key extracellular-signal-regulated protein kinase. The abnormal hyper-phosphorylation of tau, a gene that encodes the microtubule-associated protein tau, promotes microtubule stability and functions as a linker protein between axonal microtubules and neural plasma membrane components. Its mutation in Alzheimer’s disease (AD) has been shown to involve the ERK of the mitogen-activated protein (MAP) kinase family. Both the intracellular and regional distribution of the active forms of both MEK1/2 and ERK1/2, and the accumulation of p-MEK1/2 and p-ERK1/2 are known, the latter found in cases with stages I–III neurofibrillary degeneration (36).

We detected the differentially expressed ncRNAs across DRG and SN. Further classification was operated into three lncRNA categories, namely pseudogenes, lincRNAs, and asRNAs. LincRNA sequences categorized into different types on the basis of their genomic location with respect to protein-coding genes have been explored in very few studies (9, 22). The lincRNAs, detected by using both mouse and rat references, were manually curated to verify whether they could reveal potential candidate therapeutic targets for NP. In particular, potential association with the neighboring DEG targets was selectively investigated, i.e., within ±3 Mbps.

Table 2 shows the list of DEGs with respect to lincRNAs in both tissues and the corresponding annotations. The directional location of putative targets with respect to the genomic location of the lincRNAs is also reported. Analysis of lincRNAs and their neighboring DEGs revealed about 26 potential protein-coding targets in SN and 4 in DRG. In particular, we enabled a search for context-related information for each candidate bioentity, and after further manual curation of these 30 targets, we found examples, such as RGS22, RGS18 (regulators of G-protein signaling), and STIM2 (stromal interaction molecule 2), among the candidate DEGs in SN (indicated by *), with the potential to be targets of lincRNAs.

In particular, STIM2 plays an important role in ischemia-induced neuronal damage, and its absence in knockout mice was shown to interrupt blood flow in the brain, thus decreasing the neuronal damage caused by ischemia (2). The neuroprotective influence of STIM2-deficiency following an ischemic incidence suggests that the inhibitors of the STIM2 function might have potential therapeutic effect as neuroprotective agents to treat ischemic injury and other neurodegenerative disorders showing altered Ca²⁺ homeostasis. This basically suggests the role of STIM2 in hippocampus-dependent spatial memory, synaptic transmission, and plasticity. Since NP is likely a result of long-term plastic changes along somatosensory pathways, we hypothesize that STIM2 might play a role in NP-related synaptic transmission.

In addition to STIM2, we also detected a few genes from the family of RGS (regulators of G-protein signaling), namely, RGS18 and RGS22, as putative targets for lincRNAs. RGS proteins share a conserved domain of ~120 amino acids and are responsible for accelerating GTPase activity on the G-protein alpha subunit. They also affect physiological regulation of G-protein-mediated cell signaling in many tissues and organs and have been considered as potential drug targets in various nervous system diseases (37, 38). A family of RGS genes expressed in spinal cord may be involved in the regulation of GPCR signaling and adaptive changes of the nervous system, accompanying insensitivity to morphine observed in NP models (39). Based on the expression of RGS proteins, novel targets for small molecule inhibitors might provide specific treatment for various pathophysiological conditions (40).

For instance, RGS18 acts as a negative regulator of the osteoclastogenesis (41) and also plays an essential role in the regulation of megakaryocyte differentiation and chemotaxis (42). Finally, RGS proteins accelerate the deactivation of G proteins to reduce the GPCS (G-protein-coupled receptors) signaling, while few have effector function and thus transmit signals. The range of functions of RGS proteins along with their dynamically regulated distribution in brain makes them putative targets for therapeutic use (43, 44). While there are no current studies that show the relevance of RGS18 and RGS22 to NP, other components of RGS proteins have been shown to contribute to NP etiology (45, 46).

The contextual analysis of all the pseudogenes detected in both the neuronal tissues was performed. These pseudogenes sequences were mapped to the protein-coding genes to identify their associated parental genes, keeping only the top unique hits for downstream analysis. We further found the differentially expressed genes in either or both of the tissues. The detected pseudogenes and their corresponding parental protein-coding genes were listed in DRG (Table S2 in Supplementary Material) and in SN (Table S3 in Supplementary Material), both at the four different time points. These lists of entities were then used to draw STRING networks (Figure 3A), in an attempt to find the protein–protein interactions involved with pseudogene targets. The rat reference was used as a knowledge base, and a quite stringent confidence level of 0.7 was chosen to control for false positives. The interactions were generated by a database of known and predicted protein–protein interactions, following the STRING structure.

Figure 3B shows the protein–protein interactions for pseudogene targets (black boxes) obtained using mouse as the reference. Dotted lines indicate associations between pseudogene and protein-coding targets, which have been detected, but for which there is no interaction. In addition, we generated the protein–protein interactions for pseudogene targets obtained using mouse and rat references (at the same confidence level of 0.7, as before) for DRG (Figure S1 in Supplementary Material) and SN (Figure S2 in Supplementary Material). The parental protein-coding targets were further manually curated to match possible transcription factors, kinases, or receptors (shown in Table S1 in Supplementary Material). Protein-coding targets were scrutinized to confirm their role in various neurodegenerative processes. Table 3 reports network-driven identifications of parental genes appearing as paths in Figures 3A,B. Heatmaps were generated for pseudogenes and their DEG targets (Figures 4A,B, for SN and DRG, respectively).

Overall, we found that the expression of pseudogenes and their respective targets are more prominent 14 days after injury in DRG, whereas in SN, they are highly expressed at 1, 4, and 7 days after injury. Tables S4 and S5 in Supplementary Material contain the list of pseudogenes and their corresponding parental protein-coding genes in DRG and SN that are differentially expressed in DRG and SN at four different time points, respectively. We found about 12 DEGs in DRG and 134 DEGs in SN to be considered as potential targets for pseudogenes. In particular, we identified a number of olfactory receptors as targets of pseudogenes, which connect according to the network path Olr322-Olr522-Olr490-Olr951-Olr1565-Olr1174 in the rat network and are mostly expressed in SN, with the exception of the path Olr322-Olr522-Olr1565-Olr1174 expressed in DRG as well (Figure S1 in Supplementary Material). In general, these receptors have been reported to be important components for the early diagnosis of neurodegenerative diseases (56–63).

We then identified HSPD1 (Table 3, second row), a heat shock 60-kDA protein1, which encodes a member of the chaperonin family. Studies have shown increased frequency of the pathogenic variant of the HSPD1 single-nucleotide variant in a subgroup of sudden infant death syndrome (SIDS) patients (64). Yet, another study showed that the mutations in HSPD1, the gene encoding Hsp60, are related to two human inherited diseases of the nervous system, spastic paraplegia and MitCHAP60 disease. The study shows the significance of Hsp60 chaperone in mitochondrial function and its relation to the formation of the respiratory chain complex in neuronal tissues (47). We also detected HNRNPA1 in mouse network (Table 3, third row). This is a heterogeneous nuclear ribonucleoprotein A1. The levels of HNRNPA1 were found dysregulated in patients with AD, and the associated genotype is likely a risk factor for frontotemporal lobar degeneration (FTLD) among the male populations (49). This network also includes NAP1l1, which plays a significant role in the cortical neuronal differentiation. NAP1 is the nucleosome assembly protein 1, which is normally used for in vitro nucleosome assembly reactions under physiological ionic condition (65). Other studies have shown that NAP1 regulates actin nucleation by forming a complex with WAVE regulatory protein and is selectively expressed in the developing cortex (66). The cytoskeletal rearrangements in the cortical plate play an important role in the cortical neuronal differentiation (51). Yet, another study reveals that the depletion of NAP1 could result in reduced branching of motor axons in the neuronal development (67). Thus, spotting out NAP1 is a significant finding in this study, due to its potential involvement in neurodegeneration.

Naturally enough, we cannot claim any direct association between target genes and both physiological and pathological conditions highlighted in the studies reported in Table 3. Nevertheless, we stress the following two points. (a) The proteins assigned to network paths may or may not refer to significance levels under differentiated models and experimental conditions, but those identified here are indeed obtained from mapped DEGs. (b) A network context is ideal to identify hotspots (motifs, paths, etc.) whose relevance comes from superior robustness of interactions compared to simple correlation observed at gene expression level, which in turn makes inference more reliable.