A Time-Course Study of the Expression Level of Synaptic Plasticity-Associated Genes in Un-Lesioned Spinal Cord and Brain Areas in a Rat Model of Spinal Cord Injury: A Bioinformatic Approach

“Neuroplasticity” is often evoked to explain adaptation and compensation after acute lesions of the Central Nervous System (CNS). In this study, we investigated the modification of 80 genes involved in synaptic plasticity at different times (24 h, 8 and 45 days) from the traumatic spinal cord injury (SCI), adopting a bioinformatic analysis. mRNA expression levels were analyzed in the motor cortex, basal ganglia, cerebellum and in the spinal segments rostral and caudal to the lesion. The main results are: (i) a different gene expression regulation is observed in the Spinal Cord (SC) segments rostral and caudal to the lesion; (ii) long lasting changes in the SC includes the extracellular matrix (ECM) enzymes Timp1, transcription regulators (Egr, Nr4a1), second messenger associated proteins (Gna1, Ywhaq); (iii) long-lasting changes in the Motor Cortex includes transcription regulators (Cebpd), neurotransmitters/neuromodulators and receptors (Cnr1, Gria1, Nos1), growth factors and related receptors (Igf1, Ntf3, Ntrk2), second messenger associated proteins (Mapk1); long lasting changes in Basal Ganglia and Cerebellum include ECM protein (Reln), growth factors (Ngf, Bdnf), transcription regulators (Egr, Cebpd), neurotransmitter receptors (Grin2c). These data suggest the molecular mapping as a useful tool to investigate the brain and SC reorganization after SCI.


Introduction
An acute lesion of the Spinal Cord (SC) is followed by an ensemble of reactive phenomena, embracing the molecular, structural, electrophysiological and functional levels, which involve not only the lesion site, but all the Central Nervous System (CNS) areas anatomically related and even unrelated to the lesion. The molecular, anatomical and functional changes extend over a long time-span, from the trauma to the consolidated final outcome, and might be responsible for partial or complete functional recovery, but also for aberrant phenomena, such as motor patterns supporting spasticity and sensory patterns supporting chronic neuropathic pain [1]. For example, electrophysiological studies have suggested that neuronal circuits below the SC level of injury, thus deprived of supra-spinal inputs, undergo progressive and long-lasting remodeling [2], being this a possible source of neuronal dysfunctions [1,3].
On the other hand, while anatomical and functional restoration at the lesion site is strongly limited by the poor repair capability of the CNS and by the anatomical disruption that follows traumatic lesions, the overall recovery of a complex function, such as locomotion, involves the recruitment of alternative or dormant pathways and learning processes [4]. In fact, an alteration of the resting-state of sensorimotor network and an increased connectivity between motor components has been demonstrated in patients with Spinal Cord Injury (SCI) in several task-evoked functional magnetic resonance imaging (fMRI) studies [5], and neuroplasticity is currently a focus for rehabilitation approaches [6].
"Neuroplasticity" includes multiple levels of organization, from the molecular level in spinal and supraspinal circuits, to behavioral and functional performance [7]. Even though neuroplasticity is regarded as the biological substrate of rehabilitative motor training, and eliciting it represent one of the most widely used approaches to promote moderate recovery following injuries of the CNS, few preclinical data useful for translational purposes are available to support and provide experimental evidence on undergoing molecular and structural events [8].
One of the early events following SCI is a change in gene transcription in the CNS, that supports functional adaptation and compensation to the lesion, from early-activated to genes supporting long-lasting changes in neuronal function. The study of this complexity, that includes dozens of genes, their temporal expression profile in respect to the lesion and several CNS areas also far from the lesion, can be supported by "omics" technologies and bioinformatic data analysis and integration. However, most of the published data are based on available databases [9][10][11][12] or focused on SC tissue sampled at the lesion site at early time points [13][14][15], thus not considering the complex regional and temporal pattern of neural plasticity.
In attempt to fill this gap and highlight the molecular substrate of neuroplasticity in SCI, in this study we investigated changes in the expression level of genes associated with synaptic plasticity at 24 h, 8 and 45 days after contusive SCI in rats Motor Cortex (CTX-M), Basal Ganglia (BG), Cerebellum (CB) and un-injured SC, rostral and caudal to the lesion level. We adopted a data-driven approach starting from own experimental data obtained by animals fully characterized in term of clinical state, locomotor behavior and anatomy of the lesion, and a bioinformatic analysis set on stringent statistical criteria. The investigated genes include Immediate-Early (IEGs) and Late Response Genes, Long Term Potentiation (LTP) and Long Term Depression (LTD) associated genes, genes encoding for Cell Adhesion Molecules, Extracellular Matrix (ECM) Molecules, CREB Cofactors, Neuronal Receptors and Postsynaptic Density Proteins (PSD).

Lesion Characterization
The lesion was characterized for the functional outcome and anatomical landmarks. Animals were monitored for body weight ( Figure 1A), Basso, Beattie and Bresnahan (BBB) score (according to [16], Figure 1B) and recovery of the spontaneous bladder function, expressed as percentage of animals showing impairment of spontaneous urination, in which 100% correspond to the total number of animals ( Figure 1C). As expected, lesioned rats have a lower body weight gain (two ways ANOVA, p < 0.0001). Bladder function recovers between 7 and 10 days post lesion (DPL) ( Figure 1C). Locomotion, as evaluated by the BBB score having 21 as maximum score in all animals before surgery, partially recovers, thus stabilizing at 14 days. Results are presented as mean ± SEM; statistical analysis: two-way ANOVA, ** p < 0.01; **** p < 0.0001. Gait analysis was also performed by a computerized video-tracking apparatus, and the most significant parameters are presented in Figure 1D,E. Overall, lesioned rats show severe alteration in the gait pattern. For example, the step sequence regularity index, that is calculated as % index for the degree of interlimb coordination during gait (number of normal step sequence patterns (NSSP), multiplied by the number of paws and divided by the number of paw placements; RI = 100% × (NSSP × 4) / no. of paw placements) is severely altered ( Figure 1D), as also the duty cycle for both, hind and front paws ( Figure 1E,F). The duty cycle (%) is the ratio of stand time to step cycle (duty cycle = stand time/step cycle). Stand time (s) is calculated by the duration of contact of each paw with the walkway.
Step cycle (s) is calculated by the duration of two consecutive initial contacts of each paw (step cycle = stand time + swing time). These parameters show the different temporal dimensions (duty cycle), and interlimb coordination (RI) in sham and lesioned animals.
We then characterized the anatomy of the lesion at 8 and 45 DPL, and results are presented in Figure 2. Panel A indicates the site of the lesion. The SCs were serially sectioned according to the coronal plane from level T6 to level L3; alternate sections were stained with toluidine blue and HE, and the area of lesion measured at defined space intervals along the contused area. The resulting 3D reconstruction at 45 DPL is presented in Figure 2B, were the histograms express the lesioned area as % value of the total area. The schema in panel C reports the anatomical distribution of ascending (left side) and descending (right side) pathways, and the low-power micrographs in D illustrate the coronal section along the lesion extension. To better appreciate the dorsoventral extension of the lesion, SC was also horizontally sectioned at 8 and 45 DPL. The resulting 3D reconstruction is presented in Figure 2E,F, respectively. The overall analysis indicates that the lesion determines a severe disruption of the dorsal part of the SC over a rostro-caudal extension 3 mm ( Figure 2B, histogram) and a larger extension of the lesion caudally to the impact epicentrum. The ventral structures, although invested by the contusion, are less extensively destroyed and the histological analysis reveals, as expected, an intense infiltration of inflammatory cells, areas of hemorrhages and necrosis at 8 DPL ( Figure 2G,H). Phagocytic cells (gitter cells), having spherical shape with bubbly margin, reduced nucleus, lipid-laden, are also evident (arrows in H). At 45 DPL, cavitation is observed in communication with the central canal ( Figure 2I) and surrounded by an ependymal-like layer ( Figure 2J). Inflammatory and gitter cells are also present ( Figure 2K).

mRNA Expression Level of Synaptic Plasticity Associated Genes
The mRNA expression analysis has been carried out using the Qiagen RT 2 -qPCR array for the expression of genes involved in synaptic plasticity according the MIQE Guidelines [17], having reference genes automatically detected by the analysis software. The list of investigated genes is reported in Supplementary Table S1. For biological averaging and variance reduction, samples from each group were pooled for microarray experiments, since pooling dramatically improves accuracy for very small designs [18][19][20]. Five animals for each time-point and 5 control, unlesioned rats, were included in the study. Regulated genes were selected by relative Fold Change (FC) ≥ 2. To analyze the temporal evolution of gene expression, data from lesioned animals were normalized vs intact at all the time points (1,8 and 45 DPL).
We first analyzed the SC rostral and caudal to the lesion, and a distal cervical segment. Results are presented in Figure 3, were each graph reports the FC at different time-points (DPL), with their numeric values in tables. In the rostral segment, seven genes (Cebpd, Fos, Junb, Mmp9, Ngf, Timp1, Tnf ) resulted upregulated at 8 DPL, and only Timp1 remains upregulated at 45 DPL. Seven genes (Camkg2, Egr3, Gnai1, Grm3, Grm8, Igf1, Nr4a1) were downregulated at 1 DPL, and Nr4a1 and Egr3 remain downregulated at 45 DPL. In the portion of the SC caudal to the lesion, six genes resulted upregulated (Cebpb, Cebpd, Junb, Mmp9, Timp1, Tnf ) with only Timp1 overexpressed at all considered times points. The downregulated genes resulted differentially distributed across the considered times, with Egr3, Grm8 and Igf1 emerged at 1 DPL; Arc, Egr4, Ntf3 at 8 DPL; Gnai1 at 45 DPL and only Nr4a1 downregulated at all the considered DPLs.
We finally analysed the full gene set for synaptic plasticity in the BG and CB ( Figure 5). BG showed three upregulated (Ephb2, Grm2, Reln) and two downregulated (Bdnf, Egfr2) genes ( Figure 5A); while the CB analysis revealed that Ngf gene is the only one upregulated at 45 DPL and four genes were downregulated (Cebpb, Egr3, Fos, Grin2c) ( Figure 5B).

Pathway Enrichment Analysis
We then adopted a bioinformatic approach to identify gene clustering in the different areas and at the different time points. Time-series Gene Set Enrichment Analysis were performed by comparing the Pearson correlation between post-injury genes expression profiles, from the different CNS areas and three ideal expression profiles: an early profile peaking at 1 DPL, an early profile peaking at 8 DPL and a late profile peaking at 45 DPL. These time points correspond to acute, sub-acute and chronic injury phases. The complete list of pathways resulting from GSEA analysis are available in Supplementary Materials (Supplementary Tables S2-S19).
Significant results from the SC segments (False Discovery Rate (FDR) ≤0.1) are summarized in Figure 6, and relative data are presented in Tables 1 and 2 (corresponding pathways  Enrichment Table 1; SC caudal early 1-8 DPL profiles, Table 2); pathways with FDR ≤ 0.05 are highlighted. Temporal profiles are labeled as Early Profile (peak at 1 DPL) and Late Profile (peak at 45 DPL); 8 DPL pathway enrichments (FDR ≤ 0.05) are individually marked (* 8 DPL). Pathways (circles) are connected by edges (lines) based on their similarity score and grouped into labeled clusters using AutoAnnotate application v.1.3.4. Red dots denote positively associated pathways (enriched respect to control, NES ≥ 1.5) and blue dots denote negatively associated pathways (decreased respect to control, NES ≤ −1.5). Edge width denotes overlap between pathways.  The SC segment rostral to the lesion showed an early (1 DPL peak) upregulation of inflammation and differentiation pathways: interleukin-4 and interleukin-13 signaling   Table 3; BG late 45 DPL profile, Table 4; CB early 8 DPL, and late 45 DPL profiles Table 5); pathways with FDR ≤ 0.05 are highlighted. Temporal profiles are labeled as Early Profile (peak at 1 DPL) and Late Profile (peak at 45 DPL); 8 DPL pathway enrichments are individually marked. Pathways (circles) are connected by edges (lines) based on their similarity score and grouped into labeled clusters using AutoAnnotate application v.1.3.4. Red dots denote positively associated pathways (enriched respect to control, NES ≥ 1.5) and blue dots denote negatively associated pathways (decreased respect to control, NES ≤ −1.5). Edge width denotes overlap between pathways.

Discussion
The study of the spatial and temporal pattern of gene expression regulation in experimental animals is a key step in understanding the molecular bases of neuroplasticity associated to SCI, helpful to understand maladaptive plasticity supporting neuropathic pain and spasticity, and to map spontaneous recovery and neurorehabilitation interventions in a translational perspective [21,22].
In this proof-of-concept study we investigated the expression regulation of genes involved in synaptic plasticity also exploiting a data-driven approach for the analysis. We included in the study brain and SC areas outside the lesioned segment, i.e., in the CTX-M, BG, CB and SC segments rostral (T4−T7) and caudal (T12−L3) to the lesion (T9). Tissue samples were collected at 3 different time points, corresponding to acute injury phase (1 DPL), in which vascular damage, ionic imbalance, neurotransmitter accumulation (excitotoxicity), free radical formation, calcium influx, lipid peroxidation, inflammation, edema and necrotic cell death take place [23,24]; (ii) sub-acute phase (8 DPL), involving apoptosis, demyelination of surviving axons, Wallerian degeneration, axonal dieback, matrix remodeling and evolution of a glial scar around the injury site [25]; (iii) chronic phase (45 DPL), with the formation of a cystic cavity, progressive axonal die-back, maturation of the glial scar and stabilization of the functional recovery [26]. The inclusion of SC segments which are not directly involved in the primary lesion, brain areas responsible for the cortical, extrapyramidal and cerebellar movement control, as well as the inclusion of a late time point, together with the focus on synaptic plasticity genes, represents the novelty of this study, covering a still void segment in literature.
We first carefully characterized the animal model from both a functional and an anatomical point of view. Since SCI induces a long-lasting disruption of the estrous cycle [27], we used female rats in which the estrous stages were not determined. The full characterization of the lesion model is, in fact, a prerequisite for result robustness and reproducibility in all experiments on SCI [28]. In particular, it is important to determine if the lesion model induces a full anatomical disconnection between the rostral and the caudal segments of SC at the site of injury. In the contusion model used in this study, the lesion induced a severe disconnection, which preserves only the external part of the white matter in the ventral horn, where ascending spinothalamic and spinocerebellar, and part of the descending corticospinal and vestibulospinal descending tract pass [29]. This could explain the partial, spontaneous locomotor recovery observed in our rats.
We then mapped the expression regulation of synaptic plasticity-associated genes. Since gene expression regulation is a key step in neural plasticity, especially regarding long-lasting changes, we used a pathway array approach including almost 90 genes. The most relevant results, as summarized in Figure 8, are the following: (i) a different gene expression regulation is observed in the SC segments rostral and caudal to the lesion; (ii) CTX-M is the brain area showing more gene expression changes; (iii) long lasting changes in the SC include the ECM enzymes Timp1 [30], transcription regulators (Egr, Nr4a1), second messenger associated proteins (Gna1, Ywhaq); (iv) long-lasting changes in CTX-M includes transcription regulators (Cebpd), neurotransmitters/neuromodulators and receptors (Cnr1, Gria1, Nos1), growth factors and related receptors (Igf1, Ntf3, Ntrk2), second messenger associated proteins (Mapk1); long lasting changes in BG and CB include ECM protein (Reln), growth factors (Ngf, Bdnf ), transcription regulators (Egr, Cebpd) and neurotransmitter receptors (Grin2c). The pathway enrichment analysis of data derived from the SC suggests that the early profile (1 and 8 DPL) involves differentiation-, interleukin mediated-signaling paths and glutamate receptor-related path, while the late profile (45 DPL) is characterized by regulation of RNA transcription and synaptic cell signaling, with a silencing of inflammationrelated pathways.
The late time point included in this study reflect the ongoing molecular events supporting the synaptic function. In fact, while input deprivation induces an early decline of activity in the brain areas, it is followed by a reorganization of the receptive fields, probably based on a re-balance between inhibitory and excitatory synaptic changes that requires the lesion stabilization [31]. Being the expression of genes related to synaptic plasticity "activity-dependent", we assume that the abrupt interruption of axonal connections [32] can trigger rapid and long-lasting gene expression changes in the neurons originating axons of the descending/ascending spinal tracts, but also trans-synaptically [33]. At 45 DPL, same genes are differentially expressed rostrally and caudally to SCI site. In particular, only Timp1 remains up-regulated, as previously described [30], and Nr4a1 down-regulated in both segments; while Egr3 and Gnai1 are downregulated in the rostral and caudal segment, respectively. The different gene members of the EGR family of transcriptional regulators are involved in CNS function, and Egr3 is required specifically for short-term memory [34]. Gnai1 belong to the family of heterotrimeric signaltransducing molecules consisting of alpha, beta and gamma subunits, and the encoded proteins are part of a complex that responds to beta-adrenergic signals. This gene also emerged as hub gene in a protein-protein interaction network based on Gene Expression Omnibus (GEO) database, including 1, 3, 7 and 14 DPL [35]. While several preclinical and clinical studies indicated substantial differences in the anatomical outcome of the SC segments rostral and caudal to the lesion [36,37], the supporting mechanisms are not clear. The secondary injury extension is not specular respect to the lesion epicentrum. In a compression model, a different energy metabolism related to vasculature dynamics has been described, suggesting that a deficit in the glycolytic pathway accelerates the caudal degeneration; while immediate rostral degeneration is exacerbated by oxidative stress [38]. The extension of the damage outside the secondary lesion area, shows a different anatomical evolution according to the rostro-caudal gradient, including a different white vs. gray matter degeneration related to the distance of the axonal lesion from the cell body. Moreover, different sensory, motor and intraspinal pathways are differentially vulnerable [39,40]. Due to the strong correlation between anatomical and functional outcomes, a more detailed study of the rostro-caudal impact of the SCI is desirable. CNS brain areas involved in the motor function (CTX-M, BG and CB) are characterized by a profound gene expression regulation at both early and late phases. In particular, two genes were up-regulated, and 24 down-regulated, in CTX-M at the early phase; thus confirming at molecular level the electrophysiological and transcranial optogenetic mapping data, describing an initial loss of motor map (early phase) and a subsequential partial recovery (late phase) [31,41,42]. The regulated gene expression include inflammatory genes and a downregulation of Mapk2, as described at 14 DPL in a whole brain analysis [43]. In the late phase, 3 genes are up-regulated and 6 down-regulated, and the pathway analysis indicates the NMDA transmission regulation as core of the regulated net. NMDA is a glutamate receptor and ion channel found in neurons, possibly mediating the increased excitability of the CTX-M. The elective activation of glutamate neurons in the primary CTX-M was related to the functional prognosis of patients and can promote functional recovery after SCI [44]. Other regulated genes not directly related to the NMDA net includes Cnr1 (Cannabinoid Receptor 1), which is involved in sensory deprivation-induced cortical plasticity by mediating long-term potentiation and depression of synapse strength and the fine-tune of excitatory/inhibitory balance [45]. Several growth factor-and growth factor receptor encoding genes are down-regulated, e.g., Igf1, Ntf3, Ntrk2. IGF-1 affects the size of cortical receptive fields and the cutaneous threshold and is implied in cortical changes due to hypoactivity [46]; Ntf3 is among the differentially expressed genes in the auditory cortex following auditory deafferentation [47]. Ntrk2 encodes for Brain Derived Neurotrophic Factor (BDNF) high affinity receptor, TrkB, being BDNF a key growth factor for cortical neuroplasticity following lesions and physical exercise [48]. These results are not surprising, considering that learning processes involving BDNF in several spinal and supra-spinal circuits, induced by training and experience, can promote recovery after SCI [4]. In addition, human studies have described motor cortical maps reorganization after SCI [49][50][51]. These studies have also provided the background and the appropriate endpoint to evaluate efficacy of neuro-rehabilitation and neurostimulation protocols [52], evoking "neuroplasticity" as possible mechanism [7]. The gene expression regulation that we observed in our experiment possibly reflects the deep reorganization of motor pathways that starts immediately after injury [53,54].
Notably, BDNF encoding gene is upregulated at 45 DPL in the BG, while NGF encoding gene is upregulated in the CB. BDNF is a major player in the BG biology and function, also mediating various synaptic reorganization processes [55]. Moreover, BG and all extrapyramidal pathways undergo extensive reorganization changes, the magnitude of which predicts the functional recovery after SCI [56]. NGF endogenous levels increases in the CB in experimental conditions stimulating neuroplasticity (environmental enrichment) also after lesions [57], by modulating glutamatergic transmission [58].
In the context of transcriptomic mapping for the SCI, bulk RNA-seq and single cell RNA-seq (scRNA-seq) are probably the most powerful tools to build a complete picture of the complex cellular microenvironment dynamics. When bulk RNA-seq is used in a rat model of SCI analyzed at early (1, 6 DPL) and late (28 DPL) phases, the pathway enrichment analysis revealed a prevalence of immunoresponse, but the respective contribution of the resident and infiltrating immune system cells to the lesion progression can't be established [59]. More information can be obtained by scRNA-seq, but most of these studies have been performed on mouse models, which show huge differences with rats, especially in cavitation and scar formation, and is less similar to human SCI progression [60]. ScRNA-seq in mouse SCI acute (1, 3 DPL) and sub-acute (7 DPL) indicated that the main response is mediated by resident and infiltrating immune system and OPC activation [61]. This result is in line with our observations of the activation of ECM reorganization genes in the context of the synaptic plasticity pathway. Interestingly one of the main genes involved in recovery, Igf1, identified in mouse SCI scRNA-seq, resulted downregulated in the acute phase in our rat model [62].
Overall, our data suggest that a molecular mapping could be a useful tool to investigate the brain and SC reorganization after SCI, considering the need to identify reliable and simple tools to evaluate efficacy of potential therapeutic approaches in preclinical models of SCI [21,22]. "Omics" technologies might significantly impact as screening approach for target identification. However, we consider that reliability of these bioinformatic approaches should be based on own data, or on data derived from databases sharing the full anatomical and functional characterization of the model. All animals were housed in pairs in plastic cages with standard bedding and diet ad libitum. One week before surgery all animals were handled and accustomed to bladder manipulation. Five animals for each experimental time were included in each experiment.

Animals, Surgery and Care
Animals underwent a contusive spinal lesion at the vertebral thoracic level (T9). Briefly, rats were pre-medicated with enrofloxacin and tramadol (4 mg/kg, s.c.), then anesthetized with isoflurane (1-3%) in O 2 . Rats were fixed in the stereotaxic table, and an incision with a length of 4 cm was made in the skin of the back. After cutting, the muscles were dissected layer by layer to fully expose the processus spinosus of T9-T11. The processus spinosus and lamina of T9 was removed by clamp to expose the spinal canal and spinal dura. Contusive lesion of the spinal cord was obtained with Impact One impactor (Leica BioSystems, Wetzlar, Germany) using a 2 mm tip with a force of 1 N (0.75 m/s) and 0 s of stance time; the depth of impact was 2 mm to reach ventral horns of gray matter. After performing SCI lidocaine (20 mg/mL) was administered topically. Back muscles were sutured, and the skin incision closed with wound clips. Upon completion of the surgery, animals received tramadol (4 mg/kg, s.c.) for 7 days as an analgesic and enrofloxacin (4 mg/kg, s.c.) for 7 days to prevent infection. Sham animals, receiving surgical treatment, laminectomy and pharmacological treatments without SCI, were used as control. Bladder were manually expressed twice a day until automatic voiding returned spontaneously. Animals were housed in single cage for the first week after surgery.
Evaluation of rat wellness was performed by body weight monitoring and a clinical score [63], evaluated daily for the first two weeks, then once a week until the day of sacrifice. In case of infection of lower urinary tract, animals were treated twice a day with enrofloxacin (4 mg/kg, s.c.) for three days. At the time of sacrifice, the tissues of interest (cerebral CTX-M, SCI tract T8−T11, SC segments caudal to the lesion T12−L3 and rostral to the lesion T7−T4) were dissected, and immediately snap frozen and stored at -80 • C till used.

BBB Score, Locomotion and Gait Analysis
For evaluation of hind limb functional locomotor loss, BBB score [16] was performed three DPL, and lesioned animals with a score greater than 1 were discarded from the analysis (4 animals in Lesion 8 days group and 1 animal in Lesion 45 days group were excluded from analysis). The BBB is a semiquantitative scale based on locomotor response of rats, that can take on values ranging from zero (no observable movement of the hindlimbs) to 21 (normal locomotion). BBB score was repeated once a week after surgery in both lesion, and sham animal groups to assess spontaneous motor recovery.
Gait analysis was performed with CatWalk XT (Noldus, Wageningen, The Netherlands) automatized system. Animals were trained before surgery to walk repeatedly along the platform, then tested two days before SCI and once a week after lesion. All animals underwent 3 compliant runs, as defined by instrument parameters (run duration from 0.5 to 8s), for each time point, and means of all parameters were calculated by CatWalk software. Gait analysis was performed on 11 different parameters, divided in 4 categories: Spatial Parameters (Print Area, Max Contact Area, Base of Support); Kinetic Parameters (Stand Time, Swing Time, Swing Speed, Single Stance); Comparative Parameters (Stride Length, Step Cycle); Coordination Parameters (Duty Cycle, Step Sequence Regularity Index). All parameters were analyzed for both hind paws, and front paws (represented as mean of right paw and left paw).

Histology
At the day of sacrifice rats were perfused with 4% paraformaldehyde and picric acid saturated aqueous solution in 0.1 M Sörensen buffer pH 7, then spinal cord tissue was dissected and post-fixed for 24 h, then washed with Sucrose 5% O/N. 14 µm thick sagittal and coronal sections were then prepared (Leica CM1950) and processed for histochemistry staining. Toluidine blue and Hematoxilin/Eosin (HE) were performed for evaluation of lesion area and inflammatory infiltrate. To define the area of lesion, sections were captured with Nikon Microphot-FXA equipped with a CCD camera Nikon DXM1200F (Nikon) and then measured with Photoshop (v. CS6; Adobe); all images were capture with a 4× magnification of the objective and reconstructed with Photoshop's photomerge function. Lesion area was then determined for each reconstructed section with ImageJ software (NIH) as the number of pixels occupying the lesion site, and to obtain the ratio between lesion and healthy tissue, total section areas were measured for each section. 3D reconstruction was obtained aligning different levels of the same SC (sampling step 210 µm).

Total RNA Isolation, Reverse Transcription and PCR-Arrays
Right cerebral CTX-M, CB, BG, and SC segments (caudal and rostral to the lesion, and cervical) were homogenized, and total RNA isolation was performed using RNeasy Microarray Tissue Mini Kit (Qiagen, Hilden, Germany). Total RNA was eluted in RNase Free Water and using a spectrophotometer (Nanodrop 2000, Thermo Scientific, Waltham, MA, USA), absorbance values at 260, 280 and 320 were measured. RNAs were pooled for each group, to obtain a total amount of RNA of 0.5 µg and then retrotranscribed using the RT 2 First Strand Kit (Qiagen). Synaptic Plasticity analysis was performed using Rat Synaptic Plasticity PCR Array (PARN-126ZA-Qiagen), a specific panel of 84 genes involved in Rat Synaptic Plasticity (Supplementary Table S1). Genes identified as housekeeping by the default asset of the array plate were included in the analysis due to observed gene ex-pression regulation. Real time amplification of cDNA pools was achieved with CFX96 Real Time PCR System (Biorad, Hercules, CA, USA). Each array contains internal controls for genomic DNA contamination, for the quality of the amplification, and for the variability between plates. Thermal profile of PCR reactions was performed as follow: an activation step of Taq polymerase (95 • C, 10 min) followed by 40 cycles of denaturation (95 • C, 15 s) and annealing/extension (60 • C for 1 min and 30 s). At the end of the amplification cycles the dissociation curve was obtained by following a procedure consisting of first incubating samples at 95 • C for 1 min to denature the PCR-amplified products, then ramping temperature down to 65 • C and finally increasing temperature from 65 to 95 • C at the rate of 0.5 • C/s, continuously collecting fluorescence intensity over the temperature ramp. Analysis of genes expression of the array was performed using the RT 2 Profiler PCR Array Data Analysis software v. 3.5 (SABiosciences, Qiagen, Frederick, MA, USA), and gene expression was normalized on reference genes suggested by the software (Rplp1 for SC samples and Hprt1 for brain samples).

Functional Pathway and Gene Set Enrichment Analysis
Time-series Gene Set Enrichment Analysis were performed by comparing the Pearson correlation between post-injury genes expression profiles. The complete gene lists for all comparisons are given in Supplementary Table S1. The resulting normalized gene sets expression values were used as input for GSEA desktop application version 4.1.0 (build 27) and tested for enrichment in selected molecular signatures database (MsigDB) collections: H collection: hallmark gene sets v7.4, C2 collection: CP: REACTOME database v. 7.4 and C5 collection: GO: BP database v. 7.4 [64,65]. The analysis was conducted with gene set parameters adjusted to a minimum size of 5 and a maximum size of 88 genes, to correct for the number of genes being tested. Any signature with absolute NES ≥ 1.5 and FDR adjusted p-value ≤ 0.1 were reported and used as input for GSEA Leading Edge Analysis (Supplementary Figures S13-S20). Mapping of pathway terms to pathway IDs was obtained via R v. 4.04 Bioconductor v3.12 GO.db v. 3.12.1 package and additional custom scripting [66,67]. Enrichment map visualizations (Figures 6 and 7). were obtained through Cytoscape desktop application v. 3.8.2 and its 'EnrichmentMap Pipeline Collection' v. 1.1.0 application, with default parameters [68].

Statistical Analysis
Five animals per group were used for the experiments. Each data set was analyzed using two-way ANOVA, followed by post-hoc tests for exposure and genotype. Post tests were applied only for statistically significant two-way ANOVA interactions (p < 0.05). Student's t-test was used for two-groups comparison. A probability level of p < 0.05 was considered to be statistically significant.