Mapping Motor Neuron Vulnerability in the Neuraxis of Male SOD1G93A Mice Reveals Widespread Loss of Androgen Receptor Occurring Early in Spinal Motor Neurons

Sex steroid hormones have been implicated as disease modifiers in the neurodegenerative disorder amyotrophic lateral sclerosis (ALS). Androgens, signalling via the androgen receptor (AR), predominate in males, and have widespread actions in the periphery and the central nervous system (CNS). AR translocates to the cell nucleus when activated upon binding androgens, whereby it regulates transcription of target genes via the classical genomic signalling pathway. We previously reported that AR protein is decreased in the lumbar spinal cord tissue of symptomatic male SOD1G93A mice. Here, we further explored the changes in AR within motor neurons (MN) of the CNS, assessing their nuclear AR content and propensity to degenerate by endstage disease in male SOD1G93A mice. We observed that almost all motor neuron populations had undergone significant loss in nuclear AR in SOD1G93A mice. Interestingly, loss of nuclear AR was evident in lumbar spinal MNs as early as the pre-symptomatic age of 60 days. Several MN populations with high AR content were identified which did not degenerate in SOD1G93A mice. These included the brainstem ambiguus and vagus nuclei, and the sexually dimorphic spinal MNs: cremaster, dorsolateral nucleus (DLN) and spinal nucleus of bulbocavernosus (SNB). In conclusion, we demonstrate that AR loss directly associates with MN vulnerability and disease progression in the SOD1G93A mouse model of ALS.

INTRODUCTION
Amyotrophic lateral sclerosis is a progressive and fatal neurodegenerative disorder characterised by the loss of motor neurons (MNs) (1). ALS is defined by the loss of both upper and lower motor neurons. MNs are found throughout the length of the neuraxis with axonal processes spanning the peripheral nervous system to distant effector muscles. Glutamatergic upper motor neurons (UMN), referred to as corticospinal motor neurons (CSMNs), have cell bodies located in layer V of motor cortex. In humans, these large projection neurons are known as Betz cells, forming monosynaptic connections directly on lower motor neurons (LMN) in the brainstem and spinal cord. In mice, cortical to brainstem/spinal connectivity is indirect via interneurons (2). Cholinergic LMNs have cell bodies located in discrete brainstem clusters, also known as cranial nerve (CN) nuclei, and throughout the spinal cord ventral horn, known as spinal MNs.
Populations of somatic motor neurons show resistance to degeneration in ALS patients and mouse models. Two well established examples of this include, brainstem MNs innervating the extraocular muscles, and Onuf's nucleus, innervating the pelvic-perineal muscles controlling urinary and faecal continence and sexual function (3)(4)(5)(6). MNs are also differentially susceptible to degeneration based on their target fibre type innervation. The alpha (a)-MNs forming fast-twitch, fast-fatigable motor units are more vulnerable, compared to slow-twitch, fatigue-resistant motor units (7). The gamma (g)-MNs which innervate the intrafusal muscle fibres, are resistant in ALS (8). Understanding the differences in molecular signatures of these MNs and the mechanisms of vulnerability are of great therapeutic interest in ALS.
To date, there has been little exploration into how sex steroid hormones may influence MN vulnerability. It is known that clinical heterogeneity in ALS stems from the site of onset and subsequent patterns of spread (9,10). Sex-based differences in site of onset are now well established in clinical ALS (11). Bulbar onset occurs in one third of classical ALS patients and is characterised by loss of brainstem LMNs supplying oropharyngeal muscles. It is associated with a faster disease progression compared to limb/spinal onset (12) and is more common in females (11). Conversely, limb onset, which accounts for approximately two thirds of classical ALS cases, is more common in males where the male:female ratio can approach 2:1 in sporadic ALS (12). Respiratory onset, resulting from loss of phrenic MNs within the cervical spinal cord, occurs in 3-5% of cases and carries a poor survival prognosis with male bias (9,13). The ALS variant, flail arm syndrome, resulting from degeneration in cervical LMNs innervating the upper limbs, also has a high male:female ratio of 4:1 (13). Understanding a potential relationship between sex hormone action and MN vulnerability could be beneficial in understanding the mechanisms of sexual dimorphism inherent in ALS.
AR is a nuclear steroid hormone receptor, activated by androgen binding, to initiate the classical AR signalling pathway, resulting in transcriptional regulation of target genes. Upon binding testosterone or dihydrotestosterone (DHT), AR protein conformation is stabilised (14), and the complex dimerises (15) and translocates to the nucleus where it binds to androgen response elements on DNA of target genes. A hypothesis that loss or dysfunction in AR may play a role in ALS was first proposed in 1980 (16), however, findings supporting an association are not well established. An original exploration into AR distribution, assessing tritium-labelled DHT binding, reported high AR binding in the a-MNs of the lower thoracic and lumbar spinal cord, and in some CN brainstem MN populations, excluding those resistant in ALS (17). Later, Ogata et al. reported that AR immunoreactivity was evident in the resistant MNs of CNs III, IV and VI, as well as Onuf's nucleus, challenging the theory that AR absence conferred resistance in ALS (18). The discovery that an expansion mutation in the AR gene was responsible for LMN degeneration in spinal bulbar muscular atrophy (SBMA) (19), highlighted AR abundance in MNs and their selective vulnerability to aggregating disease proteins. The primary mechanism in SBMA is induction of nuclear AR aggregation upon binding of androgens. However, evidence also suggests a loss in AR function may contribute to disease pathogenesis (20,21). Although, this gain of function mutation shares similar pathogenic mechanisms with ALS, a potential overlap in AR dysfunction in MNs remains unclear.
Here, we provide a comprehensive analysis of AR expression in MN populations vulnerable in the endstage SOD1 G93A mouse. Importantly, we mapped AR expression throughout the motor neuraxis of male mice, observing a profound loss of AR immunoreactivity within MNs of SOD1 G93A mice, with spinal and bulbar populations highly expressing AR being preserved. Our assessment of MN populations in the SOD1 G93A mouse revealed that spinal MNs in the lumbar spinal cord showed greatest vulnerability, followed by cervical MNs. CSMNs and the majority of brainstem MN populations were spared, other than the trigeminal and facial nuclei which showed mild degeneration. Further analysis revealed that lumbar MNs in SOD1 G93A mice showed loss of nuclear AR, which was apparent at presymptomatic age, prior to any overt MN loss. To our knowledge, this is the first evidence directly associating AR loss with MN vulnerability and disease progression in the SOD1 G93A mouse model of ALS.

Animals and Tissue Collection
All animal experiments were conducted in accordance with the Australian National Health and Medical Research Council published Code of Practice and the ARRIVE Guidelines (22). Approval was granted by Florey Institute of Neuroscience and Mental Health Animal Ethics Committee (approval number 15-060-FINMH) to conduct this project. Transgenic SOD1 G93A mice (B6.Cg-Tg(SOD1* G93A )1Gur/J line, stock number 004435, Jackson Laboratory) were bred on a C57BL/6J background and group housed at the Florey Institute of Neuroscience and Mental Health Core Animal Services under standard 12 h light-dark conditions with access to standard rodent chow and water ad libitum. The clinical symptoms of this SOD1 transgenic mouse have been described (23). Briefly, mice develop clinical symptoms from P90 which are usually clearly evident by P120, and include impaired hindlimb splaying when suspended, hindlimb weakness with impaired gait and weight loss. For motor neuron counts and AR expression mapping, SOD1 G93A males were aged to clinical endstage. This was defined as advanced but incomplete paralysis of the hindlimbs or complete hemiparalysis in either hindlimb, and/or a decline in body condition resulting in a cumulative loss of 20% peak body weight. In the current study, the average age of endstage was 160 ± 12 days for SOD1 G93A males. Tissue from non-transgenic (wildtype) littermates was collected at P150. We selected this endstage time point to observe maximal cell degeneration in motor neuron populations outside of the lumbar ventral horn which is the most vulnerable population in the SOD1 G93A mouse. Additional cohorts of mice age P60 (presymptomatic) and P120 (symptomatic) were used for lumbar tissue collection. Mice were killed by administration of sodium pentobarbitone (100 mg/kg, i.p.) followed by cardiac perfusion with 0.1 M PBS then 4% paraformaldehyde. Dissected brains and spinal cords were postfixed for 2-4h, rinsed and cryoprotected in a sucrose gradient of 10, 20 and 30% sucrose in 0.1 M PB over 5 days. Tissue was frozen in isopentane cooled over dry ice and stored at -80°C prior to analysis.

Immunohistochemistry-DAB Staining for ChAT
Brains and spinal cords were embedded in Tissue-Tek O.C.T embedding media (Sakura Finetek, CA) and cryosectioned at 20 µm in a 1 in 10 series onto poly-L-lysine coated glass slides. Antigen retrieval was performed by baking slides for 2 h at 98°C in 10 mM citric acid, pH 6.0. Sections were blocked for 15 minutes in 0.5% hydrogen peroxide in PBS and 1 h at room temperature in 10% normal donkey serum in 0.3% Triton X-100 containing PBS. Slides were incubated in goat anti-ChAT primary antibody (1:200, Millipore, Cat# AB144P, RRID : AB_2079751) in 6% donkey serum:0.3% Triton X-100 containing PBS, at 4°C for 48 h. Donkey anti-Goat IgG HRPconjugate secondary antibody (1:500, Thermo Fisher Scientific, Cat# A16005, RRID : AB_2534679) was added to slides and incubated for 2 h at room temperature. DAB colorimetric reaction was performed using SignalStain ® DAB Substrate Kit (CST, Cat# 8059) according to manufacturer's instructions.

Image Acquisition and Analysis
Chromogenic-stained sections were imaged on a Leica DM LB2 microscope (Leica Microsystems, Germany). For brainstem cranial nerve counts, nuclei regions were captured using 20x objective tiled scan, and where possible 3 sections were counted and averaged to give single hemisphere counts for each mouse. Counts for CNs X and XII represent both hemispheres due to their medial location and the merging of hemispheres. A count from n=6 mice was averaged for each genotype. For spinal cord cervical (C1-8) and lumbar (L1-5) regions, ChAT-stained motor neurons were counted in both hemispheres of the ventral horn under 20x objective, and averaged across 28-33 sections to give a count per mouse. A ventral horn count from n=5 mice was averaged for each genotype. In cranial nerves V and VII and spinal cord sections the smaller g-MNs, were distinguished from larger a-MNs predominantly by their size, their more rounded cell morphology and more intense ChAT-staining. A retrospective comparison of MN sizes was performed and fit to Gaussian distribution. The optimal threshold to distinguish aand g-MNs based on size, as determined by the intersection of respective Gaussian distributions, was as follows: 300µm 2 and 234 µm 2 for WT and SOD1 G93A spinal MNs, respectively; 222µm 2 and 241 µm 2 for WT and SOD1 G93A CN V MNs, respectively; and 200µm 2 and 212 µm 2 for WT and SOD1 G93A CN VII MNs, respectively. To obtain phrenic MN counts, sections containing the phrenic motor column were identified by anatomical location, appearing as a MN cluster within an inter-medio-lateral position of the ventral horn of C3-5. All identifiable MNs sections were counted to provide an average count across 6 sections. The spinal nucleus of the bulbocavernosus (SNB), ventral medial nucleus (VMN), dorsolateral nuclus (DLN) and retrodorsolateral nucleus (RDLN) were located anatomically within L6 region of the lumbar spinal cord. All sections were counted and sum reported per mouse, before being averaged across n=5 mice per genotype. All counting was performed with investigator blinded to genotype.
For endstage mapping studies, fluorescent images were acquired on a Zeiss Axio 780 confocal microscope (Carl Zeiss AG, Germany). The location of the M1 motor cortex was guided by several anatomical landmarks including the corpus callosum, alignment of the lateral ventricles and anterior commissure and appearance of the third ventricle. For quantification, 1 µm Zstack images of the M1 cortex from 20 µm slide-mounted sections were acquired using 20x objective, at Bregma 1.10, 0.86 and 0.62. The area containing M1 layer V was traced using Ctip2 density as an indicator of cortical layers, and CSMN were identified as large, pyramidal shaped neurons, staining positive for MAP2 and nuclear Ctip2. Representative regions of interest (ROI) for counting layer V CSMNs and cortical layer II/III neurons are presented in Supplementary Material (Figures S1-S3). CSMNs were counted manually by two independent counters blinded to genotype and reported as number per unit of area. Counts were then averaged across n=5 mice per genotype. For AR quantification in layer V, 50 µm free floating sections were used, and for cranial nerve nuclei, 20 µm slide-mounted sections were used. Exposure parameters were optimised for each ROI in brain, brainstem and spinal cord to the brightest AR expressing MNs within an image stack (typically this was in WT mice as these had much higher signal). Image acquisition settings and ROI were kept consistent for data collection across all mice and genotypes. Z-stack images were acquired at 1 µm intervals and from an average of n=3 mice and manually counted on Zen 2.5 lite software programme (Carl Zeiss). For quantitative analysis of AR nuclear staining intensity relative to cytoplasm, this was conducted in P60 and P120 lumbar tissue. Single Z-plane images of a-MNs in L2-5 spinal cord were acquired on a Leica SP8 confocal microscope under 40x objective with experimenter blinded to genotype. Mean grey values were determined for the nuclear area and cytoplasmic area of individual motor neurons using Leica Application Suite X (RRID : SCR_013673) and an average of~35 neurons were analysed per mouse lumbar.

Data and Statistical Analyses
CSMNs, ChAT positive and AR positive CN MNs, and spinal MNs were analysed by unpaired Student's t-test comparing WT and SOD1 G93A mice. For M1 layer AR + nuclei; a and g MN subpopulation counts; P60 and P120 comparisons of AR staining intensity and AR positive nuclei, data were analysed by two-way ANOVA. Sidak's multiple comparison test was performed to compare WT and SOD1 G93A counterparts, where F-value indicated a significant genotype effect (P<0.05). All analyses were performed using GraphPad Prism 8.0 software (San Diego, CA, USA) and data is presented as mean ± SEM.

SOD1 G93A Mice Showed Regional Loss of AR Expression Within the Motor Cortex With No Loss in Layer V Projection Neuron Number at endstage
We identified CSMNs within M1 cortex of male mice using Ctip2 and MAP2 markers ( Figures 1A, B). Ctip2 is a nuclear marker most prominent in layer V pyramidal neurons (LVPNs), including CSMNs in the motor cortex (24). MAP2 was used to determine the overall neuronal architecture to select for CSMNs based on their large cell body and apical dendrite. We analysed several regions of motor cortex, tracing layer V to determine the CSMN density within this select region ( Figure S1). SOD1 G93A mice showed no detectable loss in large Ctip2 + /MAP2 + LVPNs, presumed to represent the CSMN population ( Figures 1C, D). Likewise, no evidence of somal atrophy (Figures 1C, E) or dendrite retraction ( Figure S2), indicative of cellular degeneration, was observed in SOD1 G93A mice at disease endstage. Within the motor cortex of P150 WT mice, AR appeared in all cortical layers (Figures 1B, F). Using Ctip2 nuclear marker alongside AR we were able to delineate layers II/III and V of the cortex ( Figure S3) which were observed to show the highest intensity AR nuclei. AR nuclei appeared markedly reduced in all layers, except II/III, in the SOD1 G93A mouse motor cortex ( Figure 1F and Figure S3). This was confirmed by AR + nuclear density analysis within layer II/II and V ( Figure 1G) revealing a significant effect of both genotype (F 1,8 = 14.11; P=0.0056) and cortical layer (F 1,8 = 54.36; P<0.0005) on detectable AR. AR nuclei were most abundant in WT layer II/ III neurons which was maintained in SOD1 G93A mice (P=0.2218). In contrast, a 67% reduction in AR detectable within layer V neurons occurred in SOD1 G93A mice compared to WT ( Figure 1G). AR was highly expressed in the nuclei of almost all LVPNs of WT mice ( Figure 1H; Table 1), with SOD1 G93A males displaying a reduction in LVPNs with detectable nuclear AR staining (94% vs. 73% for WT and SOD1 G93A mice, respectively; Figure 1I). This was paralleled by a similar reduction in the nuclear immunoreactivity of AR in SOD1 G93A LVPNs ( Figure 1J). In summary, CSMNs in male mice expressed high levels of AR and were resistant to degeneration in SOD1 G93A mice. These neurons displayed a moderate loss in AR by endstage disease.

Somatic Motor Neurons of the Brainstem Cranial Nerves III (Oculomotor), IV (Trochlear) and VI (Abducens) Displayed Low AR Expression Without Degeneration in SOD1 G93A Mice
MNs innervating the extraocular eye muscles are known to be preserved in SOD1 G93A mice (3). The oculomotor nerve or cranial nerve III (CN III) is the most rostrally located of the MN nuclei in the brainstem, located within the midbrain in the position of the superior colliculus ( Figure 2A). These MNs innervate most of the muscles controlling eye and eyelid movement. No change in ChAT + MN numbers within CN III  nuclei were observed in SOD1 G93A mice ( Figures 2B, C). A significant reduction in AR + nuclear expression was evident in MNs from SOD1 G93A mice ( Figures 2B, D) with AR detectable in 21 ± 4% of MNs, compared to 49 ± 5% in WT mice ( Table 1).
To explore AR expression among MN subpopulations, we used NeuN staining to identify a-MNs, distinguishable from g-MN which have low abundance or non-detectable NeuN. The loss in AR detected in SOD1 G93A mice occurred similarly across both MN populations within CN III ( Figure 2D). The a-MNs predominate in CN III nuclei (73% of ChAT + MNs) and AR was slightly more abundant in this subpopulation with 55% staining positive, compared to 31% of the total g-MN in WT mice. The trochlear nerve or cranial nerve IV (CN IV) is located caudolaterally to the CN III nuclei with the two cranial nuclei often grouped together ( Figure 2E). The CN IV innervates a single extraocular muscle of the eye controlling internal rotation. Likewise, no loss in CN IV MNs was observed in SOD1 G93A mice ( Figures 2F, G), although a significant reduction is nuclear AR   Figure 2I). Again, no loss in CN VI MNs occurred in SOD1 G93A mice ( Figures 2J, K) with very low detection of nuclear AR in both genotypes ( Figure 2L and Table 1). Taken  The trigeminal nerve or cranial nerve V (CN V) is the first of the larger multi-branched CNs with MN nuclei located in the brainstem pons, laterally to CN VI ( Figure 3A). CN V innervates eight muscles involved in chewing, biting and swallowing functions (25). In ChAT-stained sections, g-MNs were delineated from a-MNs primarily based on their smaller size (Table S1, Figure S4 in Supplementary Material), often exhibiting a rounded morphology and higher ChAT signal intensity. Genotype (F 1,20 = 56.81, P<0.0001) and MN subtype (F 1,20 = 287, P<0.0001) had a significant effect on ChAT + neuron numbers in CN V. SOD1 G93A mice lost~41% of CN V a-MNs ( Figures 3B, C), while g-MNs (making up~23 of CN V MNs in WT mice) were preserved ( Figure 3C). In fluorescent stained tissue, the g-MNs were also identified by a lack of NeuN (26); based on this delineation 82% of ChAT MNs in WT mice were positive for NeuN and 18% were negative for NeuN. In the SOD1 G93A mouse CN V, a reduction in nuclear AR levels was evident, with only 18 ± 2% of MNs having AR + nuclear staining, compared to 68 ± 3% of MNs in WT mice ( Figure 3D). This loss was preferentially in NeuN + a-MNs, with only 16% of the remaining population staining positive for AR compared to 75% of a-MNs in WT mice. In SOD1 G93A mice, 20% of NeuNg-MNs were AR + , which is a slight reduction compared to WT g-MNs where 37% were AR + . The facial nerve or cranial nerve VII (CN VII) is the largest CN with motor neuron nuclei located in the ventral aspect of the pons ( Figure 2I) innervating the muscles of the face controlling facial expression (27). Genotype (F 1,20 = 39.3, P<0.0001) and MN subtype (F 1,20 = 11.38, P<0.0001) had a significant effect on ChAT + neuron numbers in CN VII. In SOD1 G93A mice,~22% of CN VII a-MNs were lost by endstage disease with no loss in g-MN (which represented 15% of total CV VII MNs in WT mice; Figures 3E, F and Table S1). In fluorescent stained tissue, there was regional variation in ChAT and NeuN staining intensity across MN populations ( Figure 3E) and NeuN was also present in the nuclei of most small MNs. In the CNs, g-MNs can appear similar in size to smaller a-MN (28), making subtype delineation more difficult. In WT mice, NeuN + staining was detected in 93% of ChAT MNs with variable cytoplasmic and/or nuclear intensity. Therefore, we quantified AR staining in the total MN population. Again, a significant reduction in AR levels was evident in SOD1 G93A mouse CN VII, with only 26 ± 8% of MNs observed to have AR + nuclear staining, compared to 69 ± 9% of ChAT + neurons in WT mice (Figures 3E, G; Table 1).
Within the CN VII, clusters of lower expressing AR + MNs were evident in the ventral-lateral zone ( Figure 3E, inset a), while subnuclei in the dorsal-medial zones exhibited some more intensely stained AR + ChAT + neurons ( Figure 3E, inset b). To explore further the regional specific variation in MN vulnerability and AR expression, the CN VII was divided into 6 subnuclei, which are named relative to their anatomical location within CN VII nucleus ( Figure 3H), also reflecting topographical innervation of the facial muscles (29). We selected two CN VII subnuclei for further investigation; the lateral subnucleus (7VL) innervating the nasolabialis muscles, and the dorsal intermediate (7DI) subnucleus innervating the facial eye area ( Figures 3H, I). In SOD1 G93A mice, there was a 30% loss in MN within the 7VL subnucleus ( Figure 3J), whereas MNs in the 7DI subnucleus were preserved ( Figure 3L). A reduction in AR + MNs was evident in both subnuclei of SOD1 G93A mice, although this was more pronounced in the 7VL region, which was reduced from 72% to 14% AR + MNs in WT vs. SOD1 G93A , respectively ( Figure 3K), compared to 88% vs. 41% in the 7DI region ( Figure 3M). In summary, the CNs V and VII were vulnerable to MN degeneration in SOD1 G93A mice and showed a prominent loss in AR expression.
Visceral Motor Neurons of the Brainstem Vagus Nerve Complex, X (Vagus) and Amb (Nucleus ambiguus), Expressed High AR Levels and Were Preserved in SOD1 G93A Mice A CN nucleus rich in AR expression is the dorsal nucleus of the vagus nerve, also known as cranial nerve X (CN X), with MNs located in the medulla oblongata of the brainstem below the fourth ventricle ( Figure 4A). These MNs differ from the other CNs in that their projections form visceral efferents innervating the parasympathetic autonomic ganglia of the thorax and abdomen. The CN X MNs showed a tightly packed and rounded cell morphology with consistent low level NeuN expression ( Figure 4B). No decrease in ChAT + cell count was observed in SOD1 G93A mice ( Figure 4C). AR was highly expressed in the nuclei of CN X MNs irrespective of genotype, detectable in 86 ± 2% of MNs in WT and 90 ± 10% of MNs SOD1 G93A (Figures 4F, G). The nucleus ambiguus (Amb) is a group of motor neurons situated in the medullary reticular formation of the brainstem ( Figure 4J) and has a branchial MN component innervating the muscles of the palate, pharynx and larynx, as well as a visceral MN component innervating the heart (30). Amb MNs were preserved in SOD1 G93A mice at endstage disease ( Figures 4K, L)  staining. ****P < 0.0001 significantly different to WT male by two-way ANOVA with Sidak's multiple comparisons test comparing genotype for ChAT quantification; mean ± SEM, n = 6. *P < 0.05, **P < 0.01, ****P < 0.0001 significantly different to WT male by unpaired t-test for AR quantification and subnuclei quantification, mean ± SEM, n=3. Scale bars = µm units.

Somatic Motor Neurons of the Brainstem Cranial Nerve XII (hypoglossal) Expressed Low to Moderate AR levels and Showed Early Signs of Degeneration in SOD1 G93A Mice
Situated immediately ventral to the CN X is the hypoglossal nerve or cranial nerve XII (CN XII), which innervates the muscles of the tongue. No loss in ChAT + CN XII MNs was observed in SOD1 G93A mice ( Figures 4B, D). Further analysis performed on the morphology of the CN XII MNs revealed a 17% loss in soma volume in SOD1 G93A mice ( Figure 4E). A clear delineation between aand g-MNs was not evident in the CN XII MNs and NeuN + was detectable in 97% of ChAT + MNs. CN XII MNs expressed a low to moderate level of AR + which was decreased in SOD1 G93A mice, where 21 ± 5% of MNs had AR + nuclear staining ( Figure 4I), compared to 74 ± 5% of MNs in WT mice (Figures 4H, I; Table 1).
In conclusion, the CN XII MNs did not show any cell loss in SOD1 G93A mice, similar to earlier described somatic CN MNs innervating the eye muscles. These CN XII MNs displayed a low AR level and showed early morphological evidence of degeneration.

Motor Neurons in the Cervical Spinal Cord Showed Consistent AR Expression Levels and Increased Susceptibility to Degeneration
The cervical spinal cord contains several motor columns innervating the neck and mastoid musculature (C1-5) and the upper limbs (C5-8). It also contains the phrenic motor column, extending from C3-C5, innervating the diaphragm which we have located and shown in Figure 5A. In SOD1 G93A mice, 53% of phrenic MNs were lost by endstage disease (Figures 5B, C) and AR staining was notably reduced in the remaining MN pool. The ChAT + cervical spinal MNs showed a clear delineation between aand g-MN populations, where g-MNs had small, rounded morphology ( Figure 5D; Figure S4; Table  S1) and an absence in cytoplasmic and nuclear NeuN ( Figure 5E). Genotype (F 1,18 = 82.8, P<0.0001) and MN subtype (F 1,18 = 177, P<0.0001) had a significant effect on ChAT + neuron numbers in the cervical ventral horn. In SOD1 G93A mice, 54% of a-MNs were lost through C1-C8, while there was no change in g-MNs comparable to WT mice ( Figure 5F). A significant reduction in AR levels was observed in SOD1 G93A cervical MNs, with 46 ± 6% of remaining MNs having AR + nuclear staining (compared to 69 ± 3% in WT mice), although with markedly reduced intensity ( Figures 5E, G).
In remaining NeuN + a-MNs, 55% had low level AR nuclear staining (compared to 69% of a-MNs in WT mice) and AR was mostly absent in the NeuNg-MN population with only 6% having detectable AR signal (compared to 58% of NeuNg-MNs in WT mice). In summary, spinal MNs in the cervical region showed pronounced cell loss in SOD1 G93A mice, much greater than the neurodegeneration described in brainstem MNs. The Cr9 MNs did not undergo degeneration in SOD1 G93A mice ( Figure 6C) and AR immunoreactivity was highly conserved in these neurons (Figures 6B, C). Conversely, MNs present in the MMC and LMC throughout L1-5 were the most vulnerable MN population in SOD1 G93A mice ( Figures 6D-F) with 67% of the a-MNs lost through L1-L5 ( Figure 6F) and no loss in g-MNs (which made up 26% of the total ChAT + ventral horn MN population in WT mice, Table S1). A significant reduction in AR levels was evident in SOD1 G93A mice, with 55 ± 9% of remaining MNs retaining AR nuclear staining, albeit at markedly reduced intensity, compared to 81 ± 2% of MNs in WT mice ( Figures 6E, G). Consistent with that observed in cervical MNs, in SOD1 G93A mice the remaining NeuN + lumbar a-MNs, showed very low level AR nuclear expression, detected in 55% of a-MNs (compared to 85% of a-MNs in WT mice; Figure 6G). In the NeuNg-MN population from SOD1 G93A mice, 53% maintained detectable, but weaker nuclear AR (compared to 75% of g-MN in WT mice). The caudal section of the lumbar spinal cord running into the early sacral regions contains several sexually dimorphic MN populations which innervate the pelvic/perineal musculature in male mice. This includes the highly AR expressing spinal nucleus of the bulbocavernosus (SNB), located in the central area of L5-S1, and the dorsolateral nucleus (DLN) located, as the name suggests, in the dorsal lateral aspect of the ventral horn of L6-S1 ( Figure 6H). In SOD1 G93A mice, neither of these nuclei exhibited significant MN loss, compared to WT mice (Figures 6M, O) and maintained high AR expression ( Figures 6I, K) at endstage disease. The sparse neurons contained in the ventromedial nucleus (VMN; Figures 6J, N) and the monomorphic retrodorsolateral nucleus (RDLN; Figures 6L, P) were typical of the general lumbar MN population and showed significant reductions in ChAT + counts in SOD1 G93A mice, with 60% loss in the VMN ( Figure 6N) and 42% loss in the RDLN ( Figure 6P). In SOD1 G93A mice, the intensity of AR staining in MNs was observed to be decreased in both VMN and RDLN, although we did not quantify expression levels. The results from the lumbar spinal cord when taken together, showed that the lumbar MNs were the most vulnerable population in SOD1 G93A mice and displayed the greatest cell loss in the neuraxis. AR staining intensity was also profoundly decreased in the surviving MNs, although a low level was still detectable in half of these MNs. The lumbar spinal cord also contained several MN populations with high AR expression which did not undergo degeneration in SOD1 G93A mice.

AR Expression Was Decreased in Lumbar Motor Neurons of SOD1 G93A Male Mice From Pre-Symptomatic Age
With AR expression noticeably diminished throughout MNs of the CNS, we further explored this loss within the vulnerable MNs of the lumbar spinal cord over the disease course in male SOD1 G93A mice. P60 was chosen as an early pre-symptomatic time prior to MN loss and clinical onset. This is also a typical time for treatment initiation in SOD1 G93A mouse studies. P120 was also chosen to explore AR as clinical symptoms are established at this time but prior to the rapid deterioration toward endstage paralysis. For these studies, a-MNs were identified by ChAT and SMI-32 positive cytoplasmic staining and were selected from ventral horn regions through L2-L5 ( Figure 7A). L6 was excluded as this region contains high ARexpressing sexually dimorphic MN populations (32). Genotype (F 1,18 = 27.11, P<0.0001) and age (F 1,18 = 57.04, P<0.0001) had significant effects on a-MN numbers in the lumbar ventral horn ( Figure 7B). At presymptomatic age (P60), SOD1 G93A mice did not show any MN loss ( Figure 7B). Genotype had significant effects on nuclear AR intensity (F 1,18 = 32.8, P<0.0001; Figure 7C) and AR + nuclear scoring in a-MNs (F 1,18 = 22.89, P<0.0001; Figure 7D). At P60, AR nuclear staining intensity was significantly reduced by approximately 30% in SOD1 G93A MNs compared to WT MNs ( Figures 7A, C) with AR being detected in 60% of a-MNs (compared to 84% of a-MNs in WT mice; Figure 7D). At P120, WT mice MN count was comparable to P60 WT (P=0.1495), while SOD1 G93A mice showed a 32% loss in ChAT + MNs compared to WT ( Figure 7B). Age was not found to influence nuclear AR intensity, however, it did significantly contribute to AR + nuclear scoring (F 1,18 = 6.016, P=0.0246; Figure 7D). Further analysis did not find a difference between P60 and P120 WT mice nuclear AR (P=0.1594). In the remaining a-MNs of P120 SOD1 G93A mice, AR nuclear staining intensity was reduced by 28% the level in WT mice ( Figures 7A, C) with AR being detected in 50% of a-MNs ( Figure 7D). By contrast, the highly expressing AR MN population, Cr9, showed~3.8-fold higher AR staining intensity compared to other a-MNs, and was comparable between SOD1 G93A and WT at P120 ( Figure 7E).

DISCUSSION
Determining what influences MN vulnerability in ALS is critical for developing therapeutic interventions. MNs in the motor cortex, brainstem and spinal cord diverge across multiple physiological and biochemical properties. These include synaptic connectivity, composition of effector targets, gene expression profiles and protein homeostasis (33)(34)(35). AR is known to be differentially expressed throughout the neuraxis and the role of androgen signalling within various MN populations remains poorly understood. We previously reported that AR protein expression and androgen metabolising enzyme, 5a-reductase type 2, were reduced in the lumbar spinal cords of symptomatic SOD1 G93A male mice (36). Building on this preliminary evidence, we hypothesised that AR may influence MN vulnerability in this ALS model. In the present study, we comprehensively mapped AR expression in the MNs of the male mouse motor system, and in parallel assessed MN loss in the SOD1 G93A model of ALS. All ALS vulnerable MN populations consistently showed moderate AR nuclear expression in WT counterparts, which was robustly downregulated in endstage SOD1 G93A male mice. In the lumbar spinal cord, this downregulation in nuclear AR appeared to be present from presymptomatic age, prior to MN loss.

CSMNs in SOD1 G93A Mice and Cortical AR Expression
In sporadic and familial ALS patients, the large pyramidal Betz cells in layer V motor cortex show extensive vacuolisation of apical dendrites (37). Betz cells are synonymous with CSMNs in the mouse M1 cortex, however in mice, CSMNs are much smaller and difficult to identify (38). The clinical endstage in SOD1 G93A mice was determined by complete or hemiparalysis in the hindlimbs. At this stage, we did not identify a loss or change in the CSMNs. This contrasted with two studies in SOD1 G93A mice which reported significant losses in CSMNs. Ozdinler and colleagues reported up to 67% reduction in CSMNs in P120 SOD1 G93A , compared to WT, with CSMNs identified by retrograde fluorogold tracer (39). CSMNs were reportedly decreased by 58% in P120 SOD1 G93A mice, when identified using GFP-expression driven by the UCHL1 promotor (40). Interestingly, in both cases there were parallel losses in subcerebral projection neurons in other cortical areas including the somatosensory cortex. While using axonally transported, retrograde labelling with either virus or tracer remains the gold standard for identifying CSMN populations, these may only capture a specific vulnerable pool of the MN population. Factors such as 1) the integrity of the corticospinal tract, spinal grey matter or muscle at injection location; 2) accuracy and timing of administration; 3) confinement of retrograde label to the target area, may impact the effectiveness of these labelling techniques. Additionally, the UCHL1 gene is known to be downregulated in sporadic ALS motor cortex (41), hence, it may be lost in a disease-progressive manner, as is the case in several other neurodegenerative disorders (42,43). It is also possible that other cortical projection neurons were included in our count. In the mouse M1, corticospinal, corticothalamic, corticostriatal and corticocortical projections all emerge from layer V, with corticothalamic projection neurons showing similar size and locality to the corticospinal neurons (44). Taken together, accurate detection and quantification of CSMN loss in SOD1 G93A mice remains challenging, and no change was found in LVPNs projection neurons resembling CSMN in~P160 endstage SOD1 G93A male mice in the current study. AR signalling has been shown to regulate spine density in hippocampal pyramidal neurons (45) and attenuates dendritic atrophy in spinal MNs (46), in addition, providing neuroprotection to a range of neurons (47). Therefore, AR signalling may also play a role in supporting the health of CSMNs. In the present study we observed decreased AR staining in layer V of the motor cortex of SOD1 G93A mice, with moderately decreased nuclear AR in the larger LVPNs, representing CSMNs. With no evidence of CSMN degeneration in SOD1 G93A mice in the present study, this mouse model may not best reflect clinical ALS in which UMNs are also involved. It is possible that AR does not influence CSMNs, at least to the same extent as LMNs. In support of this, there is no sex bias in ALS presenting with predominantly UMN involvement (48), nor are UMNs involved in SBMA, a LMN disorder caused by an expansion mutation in AR (49).
We observed that AR staining was most abundant in the layer II/III neurons, with nuclear expression maintained in these neurons of SOD1 G93A mice. In rat cortex, AR was found to be present predominantly in the pyramidal neurons of layers II/III and V/VI in sensory and motor cortices (50), consistent with our findings in the mouse motor cortex. Retrograde tracer studies revealed that a large proportion of the layer V AR positive nuclei are associated with corticocortical projection neurons (50). Layer II/III neurons provide the major excitatory input to CSMN dendrites, including the thalamocortical inputs relaying cognitive and sensory information from other cortical areas (51). Cognitive impairment occurs in~50% of ALS patients with 15-20% classified as having frontotemporal dementia (FTD) (52). ALS-FTD occurs more frequently in women with increasing age (48) and cognitive impairment with executive dysfunction is also more frequent in women (53). Androgens regulate executive function (54) and are protective against cognitive decline, with multiple reports of androgen-deprivation therapy being associated with increased dementia in prostate cancer patients (55). Taken together, this evidence supports androgens and AR influencing cognitive impairments in ALS more strongly than providing neurotrophic support to UMNs.

Vulnerability of Brainstem Motor Neurons
We showed that CSMNs in SOD1 G93A do not degenerate by endstage disease. In the brainstem, we observed MN loss in several CN populations. MN loss was clearly evident in both trigeminal (CN V) and facial (CN VII) nuclei, in line with other SOD1 transgenic mouse models (3,56). We observed no loss in hypoglossal CN XII MNs, however, these MN displayed soma volume loss, an early indicator of degeneration. This data is in line with ALS patients, where the CNs V, VII and XII are all affected in early stages of disease (34). We showed that the high AR expressing MNs, innervating viscera, exhibited no degeneration and maintained their nuclear AR in endstage SOD1 G93A mice. Conversely, the MNs of the oculomotor (CN III), trochlear (CN IV) and abducens (CN VI) all expressed AR at lower levels which was detectable in less than half of the MN population. Our data provides insight into the conflicting theories on the role of AR in MNs which may have previously over generalised a relationship between AR level and MN vulnerability across LMN populations. Early hypotheses suggested that AR presence coincided with MN vulnerability and resistant CN MNs showed AR absence or low levels (16,17). Our findings are in support of this, whereby the LMNs which contained moderate AR levels were more likely to degenerate in SOD1 G93A mice, possibly being more susceptible to the loss in AR occurring throughout MNs of the neuraxis in disease. We observed this throughout the brainstem and spinal cord MNs, with low level AR populations in CNs III, IV, VI and XII showing no cell loss by endstage disease in SOD1 G93A mice. An alternative hypothesis proposed evidence that AR expression in MNs conferred resistance in ALS and this included CNs III, IV, VI and Onuf's nucleus (18). We dispute high AR expression in the oculomotor nuclei, although our findings support a theory that subpopulations of MNs expressing a high level of nuclear AR do show resistance to degeneration, as well as retaining their AR levels in endstage disease.
Androgens are protective to injured brainstem MNs. Facial nerve MNs show evidence of AR-mediated neuroprotection in axotomy models; testosterone permanently rescued~20% of MN in the postnatal hamster facial nerve axotomy model (57) and DHT transiently enhanced facial MN survival in the adult mouse facial nerve crush model (58). Testosterone improved regeneration of hypoglossal MNs following nerve crush injury in rats (59) with the tongue muscle being the primary target for testosterone-mediated neuroprotection following axotomy of the hypoglossal nerve (60). While these studies infer that androgens have neuroprotective actions against direct axonal damage in the CN nuclei, it is unknown if androgens confer the same protection against neurodegenerative MN death.

Vulnerability of Spinal Motor Neurons
Similar to cranial nerve injury models, androgens promote spinal MN recovery following sciatic nerve crush injury in rats (61), with evidence supporting androgen-mediated neurite growth and axonal recovery well documented (62). The SOD1 G93A mouse model best replicates spinal MN pathology and the dying back phenomenon, which has been likened to an axonopathy, whereby loss in synaptic connectivity at the NMJ may be the initiating site of the disease (63). AR is enriched in cells at the NMJ (64) and presents as a potential target site for androgen neuroprotection. In endstage SOD1 G93A rodents, there is an amplification in activity of remaining phrenic MNs to maintain stimulation to the diaphragm (65). Testosterone improves diaphragm neurotransmission, reducing fatigue during repetitive firing (66). The loss of AR within MNs may further impair their ability to mount effective compensatory mechanisms to maintain NMJ integrity in ALS.
Alongside protecting the NMJ and axonal health, the dendritic arbor of MNs is critical in regulating synaptic input, excitation and supporting MN health (37,67). The fast-fatigable motor units are most vulnerable in ALS, compared to slow and fatigue-resistant motor units and most cranial MNs, such as oculomotor. These vulnerable MNs tend to show a greater soma volume and increased dendritic complexity, compared to the resistant subtypes (33,68,69). Androgens provide neuroprotection against secondary dendritic atrophy induced by injury to surrounding somatic MNs (46,70). In the SOD1 G93A mouse lumbar MNs, maladaptive dendrite morphology occurs as early as postnatally, with degeneration in dendrites evident by pre-symptomatic P60 age (68), a time we observed a reduced nuclear AR. While most MNs exhibit a decrease in AR in endstage SOD1 G93A , larger spinal MNs with more extensive dendritic arborisations may be more susceptible to the impact of local androgen changes and AR downregulation.
Androgens such as testosterone and DHT bind to AR, prompting translocation of the complex into the nucleus, which also stabilises AR protein from degradation (14). In lumbar a-MNs, we report here a reduction in nuclear AR from as early as P60 prior to any cell death, and remaining diminished over disease course. The enzyme that primarily converts testosterone to the more potent DHT within androgen-responsive tissues, 5a-reductase type 2, is expressed in large pyramidal neurons in the brain and ventral horn MNs in male spinal cord (71,72). We previously reported that transcript levels of this enzyme were selectively reduced in the spinal cord of symptomatic SOD1 G93A male mice, while being conserved in prostate tissue (36), and may be responsible for localising AR loss to MNs of the ventral horn and LVPNs in the brain. Further evidence is provided from ALS patients who have markedly reduced DHT concentrations in CSF, while free testosterone levels remain unchanged (73). DHT administration to male SOD1 G93A mice improved motor function and survival (74). The loss in nuclear AR evident early in disease course prior to any MN loss, may suggest antecedent mechanisms contributing to MN vulnerability, rather than a consequence of disease course.
In conclusion, a robust decrease in AR expression levels within MNs was evident throughout the CNS of male SOD1 G93A mice. In the lumbar a-MNs, a reduction in nuclear AR staining was evident prior to MN loss and may reflect earlier reports of dysregulated local androgen biosynthesis in the SOD1 G93A mouse. Unique subsets of MNs displaying high intensity AR did not appear to undergo degeneration in SOD1 G93A mice. We did not observe a universal correlation between AR expression and vulnerability in ALS through the CNS, likely reflecting the complexity and diversity across different MN populations and the mild influences of steroid hormone signalling on MN survival. Larger, more branched, spinal a-MN are potentially more sensitive to decreased AR signalling and further exploration of the functional impacts of altered AR in these MNs is warranted.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by Florey Institute of Neuroscience and Mental Health Animal Ethics Committee.


s2fieldsofstudy: Biology Medicine

text:
Microarray and Bioinformatics Analysis of Differential Gene and lncRNA Expression during Erythropoietin Treatment of Acute Spinal Cord Injury in Rats

Purpose We performed a genome-wide analysis of long noncoding RNA (lncRNA) expression to identify novel targets for the further study of recombinant human erythropoietin (rhEPO) treatment of acute spinal cord injury (SCI) in rats. Methods Nine rats were randomly divided into 3 groups. No operation was performed in group 1. In groups 2 and 3, a laminectomy was performed at the 10th thoracic vertebra, and a contusion injury was induced by extradural application of an aneurysm clip. Group 1 rats did not receive any treatment, group 2 rats received a single intraperitoneal injection of normal saline, and group 3 rats received rhEPO. Three days after injury, spinal cord tissues were collected for RNA-Seq, microarray, differentially expressed genes (DEGs), Gene Ontology (GO) function enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, and protein-protein interaction (PPI) analyses. Results Compared with group 1, 4,446 genes were found to be differentially expressed in group 2. Furthermore, 99 lncRNAs were found to be changed in the injury group. The data indicate that 2,471 mRNAs were upregulated, and 1,975 mRNAs were downregulated in group 2 as compared with group 1. In addition, 45 of the lncRNAs were upregulated, and the other 44 lncRNAs were downregulated. The top 5 upregulated and top 5 downregulated lncRNAs that were different between group 2 and group 1 are shown. The top 5 downregulated and the top 5 upregulated lncRNAs that were different between group 3 and group 2 are shown. Conclusion RhEPO treatment alters the expression profiles of the differentially expressed lncRNAs and genes beneficial to the development of new treatments.

Introduction
Spinal cord injury (SCI) is a global health problem, and each year, there are 15 to 40 acute SCIs per million persons commonly caused by high falls, community violence, recreational activities, and traffic accidents [1]. SCI can lead to serious damage to the nervous system, including quadriplegia and paraplegia, which seriously affect quality of life [2]. The pathophysiological mechanisms of SCI are complex and involve ischemia-reperfusion leading to endothelial dysfunction and vascular permeability changes, which induce a cascade of inflammation and subsequent neuronal death and loss of neurological function [3].
Although much research has been devoted to elucidate the complex pathophysiological processes that follow SCI, no definitive treatments have been developed. Early decompression surgery can have a positive impact on outcomes. The only pharmacological treatment that is known to ameliorate neurologic dysfunction after SCI is methylprednisolone (MP). Therefore, new therapeutic strategies to promote functional recovery in SCI patients are necessary, and a better understanding of the cellular and molecular mechanisms of SCI may help in the development of new treatments [4].
Erythropoietin (EPO), also known as red blood cell (RBC) stimulating factor, is a human endogenous glycoprotein hormone that stimulates RBC production. The production of EPO is stimulated in a hypoxic environment, and EPO is used clinically for the treatment of anemia associated with renal insufficiency [5]. In a prior study, we reported that recombinant human erythropoietin (rhEPO) reduced apoptosis and inflammation and promoted myelin repair and functional recovery following compressive SCI in rats and that delayed treatment is equally effective [6]. To our knowledge no specific cellular and molecular studies have been undertaken to understand the mechanism by which rhEPO helps repair SCI.
Long noncoding RNAs (lncRNAs) [7] are RNA transcripts longer than 200 nucleotides that lack protein coding ability. Compared to well-studied protein-coding genes, the function of most lncRNAs has not been elucidated, even though a large number of genes have been identified [8,9]. However, recent studies indicated that lncRNAs are important regulatory molecules in the human genome, which exert their biological control in various ways [10,11]. lncRNAs have been associated with cell proliferation, survival, and differentiation, genomic stability, and chromatin remodeling [12][13][14]. With the development of high-throughput sequencing technology, more and more lncRNAs have been identified, and a recent study showed that lncRNA deregulation is an important factor in various nervous system pathologies and that it may play a crucial role in SCI [15].
Thus far, no studies have focused on the differential expression of lncRNAs and mRNAs in SCI treated with EPO. Thus, the purpose of this study was to examine differentially expressed lncRNAs and mRNAs by transcriptome sequencing (RNA-seq) in SCI tissues in a rat model treated with rhEPO. Data from this study may help the development of novel treatments for SCI.

Materials and Methods
2.1. Animals. The study protocols conformed to the Guide for the Care and Use of Laboratory Animals from the National Institutes of Health and were approved by the Animal Care and Use Committee of Southern Medical University. Nine adult male Sprague-Dawley rats (220-260 g) were purchased from the Animal Center of Southern Medical University. All rats were housed 3 per cage under temperature-controlled conditions, with a 12 h light/dark cycle, and had free access to tap water and food.

Experimental
Design. The 9 rats were randomly divided into 3 groups: group 1, blank control group; group 2, SCI group; group 3, rhEPO treatment group. All rats were continuously observed and fed for 3 days prior to the experiments.
No procedures were performed on the rats in group 1. They were fed and had free access to water. In group 2, a laminectomy was performed at the 10 th thoracic vertebra, and contusion injury was induced by extradural application of an aneurysm clip. The spinal cord was clamped for 30 seconds. Penicillin (1,200,000 U/kg, intramuscular) was given immediately after injury for preventing infection of the surgical incision. An intraperitoneal injection of normal saline (5 ml/kg) was given within 2 hours of the injury and repeated for the next 3 days. Group 3 rats received the same injury as group 2 rats. In group 3, a rhEPO intraperitoneal infusion (3000 U/kg) was given within 2 hours of the injury, and the same rhEPO infusion was given for the next 3 days. As in group 2, penicillin was given immediately after the injury.
Three days after SCI, spinal cord tissue was collected from each rat for the experiments described below.
2.3. Surgical Procedures. Rats were deeply anesthetized with an intraperitoneal pentobarbital injection (40 mg/kg) and were fixed in the prone position. Back hair on the surgical area was removed with electric shaver, and the area was disinfected with 3% iodophor. After locating the T10 spinous process, a 2 cm midline incision was made on the midline of the back from the T8 to T12 vertebrae. The overlying musculature was separated laterally, and the spinal cord was exposed by a complete T10 level laminectomy. Subsequently, the spinal cord was subjected to extradural compression with a temporary aneurysm clip (70g force; 65821T; Rebstock, Dürbheim, Germany) for 30 seconds to induce a crush injury. The surgical site was closed using nondegradable sutures after removing the aneurysm clip, and then, the closed skin incision was disinfected again with 3% iodophor. During the procedure, body temperature was maintained with a heat lamp.
After surgery, all rats received 2 ml of 10% glucose solution, tramadol hydrochloride (50 mg/kg) for postoperative analgesia, and penicillin (800,000 U/kg) by intramuscular injection to prevent infection of the surgical incision. The rats were returned to their cages after they completely recovered from anesthesia.
The rats were fed normally with food and water for 3 days, and manual bladder evacuation was performed 3 times a day. On the third postoperative day, the wound was recut and 1 cm of spinal cord tissue was taken from the exposed spinal cord and rapidly frozen in liquid nitrogen. After the sample was completely frozen, it was stored in an airtight container at -80°C until RNA extraction.
2.4. RNA Extraction. RNA was extracted from spinal cord tissue using the TRIzol method, according to the manufacturer's instructions.

Microarray
Analysis. The mRNA in spinal cord tissue samples was enriched with magnetic beads with the probe named oligo (dT). Subsequently, fragmentation buffer was added to break the mRNA into short fragments, and the mRNA was used as a template to synthesize cDNA using random hexamers. Double-stranded cDNA was then synthesized by the addition of buffer, dNTPs, and DNA polymerase I and RNase H. The double-stranded cDNA was then purified using AMPure XP beads.

Computational and Mathematical Methods in Medicine
Purified double-stranded cDNA was repaired, A-tailed, and ligated to the sequencing linker, and fragment size was selected using AMPure XP beads. Finally, PCR amplification was carried out, and the PCR product was purified with AMPure XP beads to obtain a final library.
After the library was constructed, preliminary quantification was performed using a Qubit 2.0 fluorometer (USA, Invitrogen), and the insert size of the library was subsequently detected using an Agilent 2100 bioanalyzer (USA, Agilent). After the insert was determined to be consistent with expectations, q-PCR was used to accurately quantify the effective concentration of the library to ensure library quality. After the quality of the library was confirmed, the high-throughput sequencing was conducted.
The raw reads, the data by sequencing, need to be filtered to eliminate low-quality reads in order to ensure the quality of the information analysis. The subsequent data obtained is recorded as total data. By removing the known ribosomal RNAs from the total data (28S rRNA, 18S rRNA, 12S rRNA, 5.8S rRNA, and 5S rRNA), high-quality, clean data was obtained. The ribosomal RNA that was removed was identified using the database of the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov). Subsequent microarray profiling was performed by the Boyue Biotechnology Company (Wuhan, China).

qRT-PCR Validation of Microarrays.
To confirm the repeatability of the microarray assays, 6 additional rats were divided into 2 groups and treated the same as the rats in group 2 (SCI group) and group 3 (SCI+rhEPO group), respectively. After total RNA was extracted from the spinal cord of the 6 rats, qRT-PCR assays were performed. Briefly, the qRT-PCR assays consisted of 2 steps, RNA reverse transcription (RT) and qPCR detection. First, the PrimeScript™ RT reagent Kit with gDNA Eraser (TAKARA) was used to synthesize cDNA according to the manufacturer's instructions after removing the genomic DNA. RT-PCR was then performed using SYBR® Premix Ex Taq™ II (TAKARA). The reaction system consists of 10 μl SYBR® Premix Ex Taq™ II, 0.4 μl PCR Forward Primer (10 μM), 0.4 μl PCR Reverse Primer (10 μM), 2 μl cDNA, and 7.2 μl ddH 2 O. The primer sequences were designed and synthesized in the laboratory by the Guangzhou cm biotechnology and are listed in Table 1. The reaction conditions were: 95°C for 10 minutes; a total of 40 cycles at 95°C for 15 seconds and 60°C for 20 seconds. Each sample tested in triplicate. Gene expression levels were normalized to glyceraldehyde 3-phosphate dehydrogenase (GAPDH) using the ΔΔCT method. Finally, Student's t-test was used to examine differences between the 2 groups, and values of p < 0:05 were considered statistically significant.

Differentially Expressed Gene (DEG) Analysis and Gene
Ontology (GO) Enrichment Analysis. DEG analysis refers to the identification of genes with significant differences in expression levels between different sample groups. The clean data was analyzed using the DESeq2 package (http://www .bioconductor.org/packages/release/bioc/html/DESeq2 .html) in the R programming language (version 3.60) because our sample had only 3 per group [16]. When the biological repeat reaches 5 to 10, a better choice is to use a nonparametric method.
DESeq2 is the most popular statistical method to analyze DEGs, and it can estimate variance-mean dependence in clean data, and test for differential expression based on a model using a negative binomial distribution. The log2fold change (log2FC) and p values of the genes were calculated [17]. A gene was considered to be a DEG when the log2FC was >1 and the p value was <0.05. The lncRNA information was then extracted from DEG analysis result based on the lncRNA annotation information provided in the reference genome annotation file.
GO enrichment analysis is used to annotate genes and gene products and to provide gene function classification labels and background knowledge of gene function and has become a common approach for sequencing data processing [17][18][19]. GO enrichment analysis can be divided into 3 parts: molecular function (MF), biological process (BP), and cell composition (CC). GO annotation information of genes can be found by searching the GO database by species and genetic information. Based on the GO annotation of a gene, all of the genes of the species can be selected as background genes, and p values can be calculated using statistical methods. As such, distribution information and the significance of the gene collection of the GO category can be obtained.
To gain further insights into the changes of biological pathways in the cells of SCI rats, Kyoto Encyclopedia of Genes and Genomes (KEGG analysis) was performed. KEEG is a database of biochemical reactions, signaling pathways, metabolic pathways, and biological processes and can be used to identify the significant pathways associated with DEGs.
The clusterprofiler package (https://bioconductor.org/ packages/release/bioc/html/ clusterProfiler.html) in the R programming language (version 3.60) allows 2 methods of analysis. In the analyses in this study, a value of p < 0:05 was considered to indicate statistical significance.
Computational and Mathematical Methods in Medicine 2.8. Protein-Protein Interaction (PPI) Network Construction. The Search Tool for the Retrieval of Interacting Genes (STRING, http://string-db.org/) is a database that aims to provide a critical assessment and integration of PPI, including direct (physical) and indirect (functional) associations [20]. In the PPI network built using the STRING online tool, each node signifies a gene, and the edges indicate interactions between nodes. The degree is defined as the number of edges linked to a given node.
Cytoscape is software that provides data integration and network visualization [21], especially in respect to the processing of databases of protein-protein interactions [22]. In the current study, the cytoHubba plug-in of Cytoscape (version 3.61) was used to screen the hub genes from the PPI network, and a node degree of ≥10 was screen as the hub genes from the PPI.

Screening and Identification of DEGs.
DEGs between samples were selected by differential multiples (log 2FC > 1 ) and a significance level of p < 0:05. Compared with group 1, 4,446 genes were found to be differentially expressed in group 2. Furthermore, 99 lncRNAs were found to be changed in the injury group. The data indicate that 2,471 mRNAs were upregulated, and 1,975 mRNAs were downregulated in group 2 as compared with group 1. In addition, 45 of the lncRNAs were upregulated, and the other 44 lncRNAs were downregulated (Table 2). A volcano plot was created that visually showed how the expressions of the lncRNAs and mRNAs changed dramatically (Figure 1).
A comparison of group 3 with group 2 identified 228 DEGs. Furthermore, 14 lncRNAs were found to be changed in group 3 as compared with group 2 ( Table 2).
Cluster analysis of DEGs was used to determine the clustering pattern of DEGs in the different groups. The top 50 genes with the largest variance in expression between the different groups, including the known lncRNAs and mRNAs, were used for cluster analysis (Figure 2).
As shown in Figure 3, the expression patterns of the selected mRNAs were consistent with the mRNA-seq results (p < 0:05 for each mRNA, Student's t test). These results indicated that the microarray were highly reliable.

Top Differentially
Expressed lncRNAs between the Groups. The top 5 upregulated and top 5 downregulated lncRNAs that were different between group 2 and group 1 are shown in Table 3. The top 5 downregulated and the top 5 upregulated lncRNAs that were different between group 3 and group 2 are shown in Table 4.

GO and KEGG Pathway Enrichment
Analyses. The gene function enrichment analysis was divided into 2 steps: gene function annotation and enrichment analysis. The 3 GO domains used to describe the gene product attributes were MF, BP, and CC (molecular function, biological process, and cell composition, respectively). GO analysis was performed for the DEGs between group 2 and group 1, and the top 15 are shown in Figure 4(a). As shown in the figure, the biological processes of the DEGs were primarily associated with positive regulation of the immune response, regulation of vesicle-mediated transport, regulation of leukocyte activation, wound healing, regulation of transmembrane transport, and other significant biological processes in SCI. In addition, the main 15 KEGG enrichment pathways were related to Epstein-Barr virus infection, focal adhesions, the calcium signaling pathway, retrograde endocannabinoid signaling, osteoclast differentiation, platelet activation, and systemic lupus erythematosus (SLE) and were highly significantly correlated with SCI ( Figure 5(a)).
Furthermore, GO and KEGG enrichment analyses were performed between group 3 and group 2 to gain further insights into the changes of biological pathways associated with rhEPO treatment of SCI. GO enrichment analysis reveals 15 significant metabolic networks, including leukocyte chemotaxis, myeloid leukocyte migration, granulocyte migration, granulocyte chemotaxis, neutrophil migration, and other biological processes that were significantly enriched in group 3 (Figure 4(b)). KEGG enrichment analysis indicated that the enriched DEGs were associated with cytokine-cytokine receptor interaction, the chemokine signaling pathway, human cytomegalovirus infection, IL-17 signaling path, malaria, and herpes simplex infection ( Figure 5(b)).

PPI Network Construction.
Based on data from the STRING database with medium confidence (data chosen had a minimum required interaction score of >0.4), a PPI network with 809 nodes and 5,081 edges was constructed between group 2 and group 1 and was visualized using Cytoscape ( Figure 6).
In the PPI network, 10 nodes were selected as hub genes using the Maximal Clique Centrality (MCC) method, which is available in the cytoHubba plug-in of Cytoscape. The hub genes were fibronectin 1 (FN1), protein tyrosine The number of this table represented the quantity of differential expression of lncRNAs and mRNAs.  Figure 1: In the volcano diagram, each point represents a gene, and the X-axis represents the logarithm of the multiple of the difference in expression of a certain gene in the two samples; the Y-axis represents the statistically significant negative logarithm of the gene expression change. The larger the absolute value of the X-axis, the greater the fold change in expression between the two samples; the larger value of Y -axis, the more significant the differential expression, and the more reliable the DEGs obtained by screening. The blue dots (fold change < −1 ) in the figure represent downregulated DEGs, the red dots (fold change > 1) represent upregulated DEGs, and the grey dots represent non-DEGs.
Similarly, a comparative study of group 3 and group 2 was conducted using the same methods. The constructed PPI network had 70 nodes and 129 edges (Figure 7). The top 10 key genes are shown in Table 6.

Discussion
In most cases, SCI is a disabling and irreversible disease that is associated with great social and economic cost to families and society [23]. SCI is a complex biological process that includes both primary and secondary damage and involves the nervous, immune, and vascular systems [24]. The initial mechanical trauma can lead to neuron necrosis and apoptosis, and the secondary damage often worsens the injury [25,26]. Active research of SCI has made slow but consistent progress with respect to developing new treatments. With respect to mRNA and lncRNA, a recent study by Jin et al. [27] identified significant DEGs at 3 days, 2 weeks, and 1 month following SCI. Based on the results of Jin's study, we choose 3 days post-SCI as the time point of our study.
Our results demonstrated that the expressions of certain mRNAs and lncRNAs were dramatically changed 3 days after SCI. Furthermore, we identified 10 key genes in the injury group (group 2) and the rhEPO treatment group (group 3) that may play an important role in early acute phase of SCI. These genes may assist in further research of SCI and the development of new treatments.  Figure 2: In the heatmap, the red color represents upregulated mRNAs or lncRNAs, and the blue color represents downregulated mRNAs or lncRNAs.

Computational and Mathematical Methods in Medicine
The genes in the injury group (group 2) with significant changes in expression profiles were associated with the positive regulation of the immune response, regulation of vesicle-mediated transport, regulation of leukocyte activation, wound healing, regulation of transmembrane transport, and other significant biological processes. Previous studies have indicated that the primary cellular responses to SCI are inflammation and an immune response, which is consistent with our GO analysis results. SCI induces the activation of immune cells and inflammatory mediators, but the benefit of targeting the immune response to treat SCI is not clear [28].
Stimulation of the proliferation of leukocytes is also a key process of SCI that occurs from the immediate phase to the chronic repair phase. Our KEGG enrichment analysis demonstrated that the genes with changes of expression were involved in Epstein-Barr virus infection, focal adhesions, the calcium signaling pathway, retrograde endocannabinoid signaling, osteoclast differentiation, platelet activation, and SLE. In a prior study, KEGG enrichment      9 Computational and Mathematical Methods in Medicine analysis revealed that the toll-like receptor signaling pathway, p53 signaling pathway, MAPK signaling pathway, and Jak-STAT signaling pathway were related to SCI [29]. Obviously, our findings are not completely consistent with that of the prior study. We postulate that the different results may be because specimens were collected at different time points. Our PPI network analysis, however, identified FN1, PTPRC, CD44, CDC20, TYROBP, AURKB, TlR2, angiotensinogen, AGT, RAC2, and MMP9 as the top 10 high-degree hub nodes, suggesting these genes may play an indispensable role in the pathophysiological processes of SCI.
Our previous studies showed that EPO reduces apoptosis and inflammation and promotes myelin repair and functional recovery following compressive SCI in rats from the perspective of bioinformatics. Our prior results were the basis for performing the current study to identify genes differentially expressed after SCI, and after treatment with rhEPO. In the rhEPO treatment group (group 3), the DEGs were significantly associated with leukocyte chemotaxis, myeloid leukocyte migration, granulocyte migration, granulocyte chemotaxis, and neutrophil migration. These findings confirm that the inflammatory response and inflammatory cell activation play an indispensable role in the repair of SCI. A prior study showed that blocking inflammation via the administration of anti-inflammatory drugs can reduce inflammation and partially restore locomotor activity following SCI [30].
Our KEGG pathway analysis indicated that the most significant pathways associated with SCI and repair were cytokine-cytokine receptor interaction, chemokine signaling, and IL-17 signaling pathways. Previous studies have reported that cytokines play an important role in central nervous system (CNS) immune system interactions, and SCI can initiate immune responses characterized by the synthesis and release of chemokines and cytokines [31,32]. In addition, an increase of IL-17 concentration can result in the increase in the size of a lesion after SCI. Study has shown that reducing the expression of IL-17 and IL-17-related inflammatory factors can protect neurons and promote recovery after SCI [33]. In the PPI network developed in this study, the top 10 high-degree hub nodes (Ccl4, Ppbp, Cxcl13, Ahsp, Ccl5, Alas2, Npy, Gng13, Ccl2, and Hba2.) were all chemokines, which are a superfamily of secreted proteins involved in immune-regulatory and inflammatory processes. Study has shown that CXCL13/CXCR5 signaling PPI network Figure 6: Interaction of protein-protein network analysis of DEGs between group 2 and group 1. Nodes represent DEGs. Lines indicate interactions between DEGs.  1  Fn1  108  6  Aurkb  61  2  Ptprc  85  7  Tlr2  60  3  Cd44  70  8  Agt  59  4  Cdc20  68  9  Rac2  58  5  Tyrobp  66  10  Mmp9  57 10 Computational and Mathematical Methods in Medicine can promote diabetes-induced tactile allodynia through the production of proinflammatory cytokines in the spinal cord of male mice [34]. Additionally, Ppbp is a platelet-derived growth factor that belongs to the CXC chemokine family. Chio et al. [35] demonstrated that a traumatic brain injury can upregulate the expression of Ppbp in peripheral blood. Alas2 encodes a protein called heme that catalyzes the first step in the heme biosynthetic pathway and appears to promote a concurrent increase of neutrophilic metamyelocytes and mature CD71 erythroid cells [36]. Furthermore, endogenous neuropeptide Y (NPY) and the activation of its associated receptors can exert long-lasting spinal inhibitory control of neuropathic pain [37]. We speculate that EPO may influence the expression of NPY to achieve behavioral NPY-induced antinociception. Taken together, the aforementioned studies and our results lead us to hypothesize that EPO plays an important role in the acute phase of SCI by regulating the inflammatory response and affecting the synthesis and release of inflammatory factors and/or chemokines. To our knowledge, this is the first study that has used bioinformatics methods to investigate EPO and the treat-ment of SCI. While we identified DEGs and lncRNAs associated with SCI, and EPO treatment, there are some shortcomings of this study. First, the small numbers of specimens may affect the reliability of the results. However, we reduced individual differences by mixing different samples from the same group. While we identified DEGs, we did not explore their functions, nor did we explore the specific mechanisms of the lncRNAs identified. In the acute phase of SCI, primary injury can lead to neuron death, demyelination in the spinal cord, and eventually, axonal dieback. In this study, we only examined tissue specimens at 1 time point. Future studies should examine and compare specimens at multiple time points after SCI.

Conclusion
In conclusion, we identified differential expression profiles of mRNAs and lncRNAs in spinal cord samples in the subacute phase following SCI. We also identified DEGs between normal spinal cord and injured spinal cord and between injured spinal cord and SCI treated with rhEPO. Critical pathways affected by rhEPO treatment of SCI were also identified. These results may offer new insights into the cellular and pathophysiological processes involved in SCI and insights for the development of new treatment methods.

Data Availability
No data were used to support this study.

Ethical Approval
The present study was approved by the ethics committee of the Zhujiang Hospital of Southern Medical University  11 Computational and Mathematical Methods in Medicine (Guangzhou, China). All applicable international, national, and/or institutional guidelines for the care and use of animals were followed.