Massive transcriptome sequencing of human spinal cord tissues provides new insights into motor neuron degeneration in ALS

ALS is a devastating and debilitating human disease characterized by the progressive death of upper and lower motor neurons. Although much effort has been made to elucidate molecular determinants underlying the onset and progression of the disorder, the causes of ALS remain largely unknown. In the present work, we have deeply sequenced whole transcriptome from spinal cord ventral horns of post-mortem ALS human donors affected by the sporadic form of the disease (which comprises ~90% of the cases but which is less investigated than the inherited form of the disease). We observe 1160 deregulated genes including 18 miRNAs and show that down regulated genes are mainly of neuronal derivation while up regulated genes have glial origin and tend to be involved in neuroinflammation or cell death. Remarkably, we find strong deregulation of SNAP25 and STX1B at both mRNA and protein levels suggesting impaired synaptic function through SNAP25 reduction as a possible cause of calcium elevation and glutamate excitotoxicity. We also note aberrant alternative splicing but not disrupted RNA editing.

abnormal calcium metabolism, altered axonal transport and activation of proteases and nucleases contribute to the pathobiology of the disease 5 .
Since motor neurons cannot be sampled during life, investigations into ALS neurodegenerative process are carried out in neuronal-like cells in culture, animal models carrying mutant transgenes or post-mortem material from patients. In the last case, the spinal cord represents the elective tissue and a pivotal source of RNA for gene expression profiling in ALS. While mixed cell populations in which astrocytes, microglia and oligodendrocytes may obscure the distinctive expression profile of motor neurons, which are the most affected in ALS, the co-culture of motor neurons and glial cells has shown an intricate interplay, suggesting that non-neuronal cells play an active role in neurodegeneration 6 . Accordingly, transcriptomic studies using high throughput sequencing approaches (RNA-Seq) to target specific tissues, such as ventral horns of the spinal cord, may provide new insights into genes and pathways involved in the degeneration of motor neurons. Indeed, RNA-Seq allows more accurate and complete expression profiling than classical microarrays, overcoming their reliance on pre-designed probes. Recently, two studies exploring transcriptome alterations in ALS by the RNAseq technology have been reported 7,8 . The first one investigated transcriptome changes in cerebellum and frontal cortex from familiar (fALS) (with mutated c9orf72 gene) and sporadic ALS (sALS) donors, showing extensive dysregulation of alternative splicing and alternative polyadenylation in both ALS groups 7 . The second one, based on RNAseq technology, investigated transcriptome changes in cervical spinal cord tissues from sALS donors showing that neuroinflammation mediated by TNF (tumor necrosis factor) is the main transcriptome abnormality in sALS 8 .
In the current work, we also report the analysis of whole transcriptome sequencing data from matched control and sALS post-mortem donors but, differently from previous studies, we focus on the detection and characterization of significant changes in ventral horns of the human lumbar spinal cord, a tissue never investigated at single nucleotide resolution. While massive sequencing of the low molecular weight RNA fraction showed 18 dysregulated miRNAs involved in TGF-beta signaling, axon survival and outgrowth, RNA-Seq data revealed altered expression of 1142 genes belonging to a restricted range of ALS deregulated molecular pathways including an impaired synaptic function. In particular, we focused our attention on two downregulated neuronal t-SNARE proteins, SNAP25 and Syntaxin-1B (STX1B), as they form with the v-SNARE protein synaptobrevin, a ternary complex mediating the neurotransmitter release at the presynaptic membrane. Notably, SNAP25 negatively regulates native voltage-gated calcium channels (VGCCs) in glutamatergic neurons 9 , and its reduced expression enhances glutamatergic neurotransmission 10 and causes neurodegeneration 11 .
Taken together, our results may represent a remarkable breakthrough to unveil sALS, supporting a primary role of t-SNARE protein expression as possible cause of intracellular calcium elevation and glutamate excitotoxicity 12 and suggesting SNAP25 and related SNARE proteins as candidate novel indicators and signatures of motor neuron death to be monitored in ALS patients.

RNA-Seq data analysis.
Our study includes ventral horns of the lumbar spinal cord from six human donors affected by sporadic ALS and five age, sex and ethnicity matched controls (Table 1). ALS and control transcriptomes were sequenced employing the Illumina technology on the NextSeq 500 platform. We generated 1.9 billion 100 bp strand-specific paired end reads with, on average, 179 million pairs of reads per sample ( Table 2). About 73% of these reads were uniquely aligned to the reference human genome (hg19 assembly) after a cleaning step to remove low quality regions ( Table 2). Using RefSeq annotations we verified that about 98% of mapped reads maintained the correct strand orientation and, on average, 46% of mapped bases corresponded to UTR and coding regions of known mRNA transcripts (Table 2).
To confirm the sporadic disease nature of our samples, we checked for the presence of known ALS mutations in aligned RNAseq data using REDItools 13 and did not find any potentially fALS-associated variant in the causative genes annotated in the ALS Online Genetic Database (ALSoD http://alsod.iop.kcl.ac.uk/) 14 . We calculated RNAseq based gene expression levels in FPKM and performed unsupervised hierarchical clustering to detect similarities in gene expression profiles between ALS and control samples. We observed a clear separation between diseased and normal groups, suggesting specific transcriptome signatures for ALS (Fig. 1A). Principal component analysis (PCA) of variant genes re-marked such separation (Fig. 1B).
Using cell type specific genes for neurons, motorneurons, astrocytes, oligodendrocytes and microglia collected from the literature 6,15 , we verified the cell type composition of our spinal cord samples ( Fig. 2 and Supplementary  Table 1). Interestingly, we found statistically significant changes in cell populations composed by neurons, motorneurons, oligodendrocytes and microglia. No substantial differences were detected in the astrocyte population (Fig. 2). While neuron, motorneuron and oligodendrocyte populations appeared depleted in ALS samples, microglial population was over-represented suggesting neuroinflammation, that is thought to be critically important in the neurodegenerative process 16 (Fig. 2).
Differentially expressed genes in ALS and controls. Differential gene expression between ALS and control samples was performed using only unique and concordant RNA-Seq reads. To improve the computational detection of differentially expressed (DE) genes, we applied two independent tools: CuffDiff2 based on FPKM metric that accounts for both uncertainty resulting from read mapping ambiguity and cross-replicate variability 17 , and DESeq2 based on counts of reads in genes 18 . The final list of DE genes was selected by the intersection between CuffDiff2 and DEseq2 lists of significant DE genes and taking into account only genes with |log2(FC)| higher than 1. On the whole, we identified 1142 deregulated genes, for some of them a specific cellular origin is known from the literature (Supplementary Table 2), with log2(FC) values highly correlated between CuffDiff2 and DEseq2 (Pearson r = 0.98) ( Fig. 3 and Supplementary Table 2). Of these, 559 (49%) were upregulated in ALS samples and mainly involved in neuroinflammation or activation of immune response (Fig. 4, Supplementary Table 3). Fifty-seven upregulated genes were also found deregulated in a previous study based on peripheral blood of ALS patients 19 and, thus, they may represent a precious source of potential biomarkers. Downregulated genes account for the 51% (583) of all DE genes and were significantly enriched (P ≪ 10 −4 hypergeometric test) in neuronal-specific genes ( Fig. 4 and Supplementary Tables 4 and 5). Consistent with past investigations, down expressed genes were involved in several biological processes of nervous system as impulse transmission (SLC1A2, CHAT, SLC12A5, HTR2C) or synapse function (SNPH, SYT4, SNAP25, STX1B) or calcium metabolism (GRIN1, GRIN2A, CACNA1G) ( Fig. 4 and Supplementary Tables 2 and 4). RNA-Seq analysis revealed also downregulation of key enzymes in cholesterol biosynthesis including HMGCR (3-hydroxy-3-methylglutaryl-CoA reductase), HMGCS1 (3-hydroxy-3-methylglutaryl-CoA synthase 1),   Within the group of down-regulated genes we found a set of 15 synaptic proteins related to the synaptic vesicle cycle (P ≪ 0.0001 hypergeometric test) that consists of docking and fusion of the vesicles, and subsequent exocytosis and neurotransmission (Table 3). Changes in synaptic transmission have been investigated in Alzheimer's disease 20 and may have a primary role in ALS for which little is known. Our results indicate deregulation of key genes for synaptic function including the two neuronal t-SNARE proteins, SNAP25 and STX1B. Recent evidence suggests novel functions for SNAP25 in controlling calcium dynamics and enhancing glutamate release in neurons 9, 10 , consistent with a role for SNAP25 in ALS pathogenesis.
SLC1A2 which encodes a glutamate transporter essential for terminating the postsynaptic action through the glutamate clearance, appeared negatively regulated in ALS. This finding was also consistent with previous investigations showing down regulation of human SLC1A2 protein in ALS spinal cord 21 .
Fourteen out of 126 genes associated with ALS (mostly by mutation or GWAS studies) from the specialized ALSoD database 14 showed differential expression in our ALS samples, 10 were down-regulated (including PRPH, UNC13A, SLC1A2, SNCG and NEFH) and 4 were up-regulated (SOD2, APOE, LOX and RNASE2).
Mitochondrial dysfunction is common in ALS pathogenesis since mitochondria have a central role in intracellular energy production, calcium homeostasis and control of apoptosis 22 . Indeed, we found several differentially expressed genes correlated with mitochondrial functions. However, the expression levels of genes transcribed by   Table 1). Statistically different cell type composition is calculated by Mann-Whitney U test (two-tails). * P < 0.05; ** P < 0.01; − Not significant. the mitochondrial genome were unchanged, indicating that mitochondrial transcription may not be impaired in ALS (Supplementary Table 6). Nonetheless, mutations in the mitochondrial genome may occur due to oxidative stress -a hallmark of motor neuron injury 22 .

Pathway analyses.
Differentially expressed genes were further analyzed through the IPA tool (Ingenuity ® Pathway Analysis, www.ingenuity.com) in order to identify the most significantly up-regulated and down-regulated canonical pathways associated to neurodegeneration in ALS.
Consistent with other transcriptome investigations in transgenic SOD1 mice and human samples, the most highly activated pathways were related to neuroinflammation and immune response ( Fig. 5 and Supplementary  Table 7). Indeed, our dataset included many up-regulated genes encoding inflammatory mediators such as chemokines and interleukins and several components of the complement system, oncostatin M signaling and STAT3 pathways (Supplementary Table 7).
Significantly down-regulated pathways, instead, included cholesterol biosynthesis, calcium signaling, glutamate receptor signaling and GABA receptor signaling ( Fig. 5 and Supplementary Table 7). In addition, several under expressed genes in ALS samples were associated with neurological and neurodegenerative human disorders such as bipolar disorder, schizophrenia, Huntington's disease or Alzheimer's disease (Supplementary Table 8). IPA analysis revealed also a significant enrichment for amyotrophic lateral sclerosis signaling, involving nine genes encoding mainly for calcium channels and glutamate receptors (Supplementary Table 8).
Our results were further corroborated by GO-PCA software implementing an unsupervised method to explore gene expression data using prior knowledge in the form of gene ontology (GO) annotations 23 . GO-PCA is able to detect specific signatures, consisting of small set of genes that are strongly correlated by their expression as well as functionally related by GO annotations 23 . We found 10 significant signatures (P ≪ 0.001) containing between 5 and 10 genes (Table 4). Five signatures, including down-regulated genes, were derived from GO terms related to cholesterol biosynthesis, glutamate secretion, motor activity or axon functionality (Table 4). Remaining five signatures, instead, were associated to inflammation and comprised up-regulated genes (Table 4).   Table 3. DE genes linked to synaptic function. List of DE genes linked to synaptic function. Per each gene we report: gene name, mean expression levels (in FPKM) in control and ALS groups (from CuffDiff2), log2(fc), P-value and corrected P-value calculated by CuffDiff2 and log2(fc), P-value and corrected P-value calculated by DESeq2.
qRT-PCR validation. We validated RNA-Seq results performing quantitative Reverse Transcription PCR (qRT-PCR) assay using RNA extracted from the ventral horns of the lumbar spinal cord in an enlarged sample group (Table 1). We included donors assayed by RNA-Seq (5 controls and 6 ALSs) and an independent group of 2 controls and 3 ALSs (Table 1). From the list of differentially expressed genes, we selected 16 candidates including genes known to be deregulated in ALS such as INA, HECW1 and SLC1A2 (representing a sort of internal control) and genes not yet described as altered in ALS but associated to significantly disrupted pathways such as the cholesterol biosynthesis, the complement system or the synaptic transmission (Table 5). All genes belonging to deregulated pathways were randomly selected (with |log 2 FC| ≥ 1), except for SNAP25 and STX1B involved in the t-SNARE complex mediating the neurotransmitter release at the presynaptic membrane.   The qRT-PCR analysis confirmed the differential expression emerged in RNA-seq analyses for all 16 genes (Fig. 6A) and the estimates of fold change in expression level were highly consistent with those from RNA-Seq (r = 0.98, P = 8.3 × 10 −11 for CuffDiff2 and r = 0.96, P = 4.2 × 10 −9 for DESeq2) (Fig. 6B).
Aberrant Alternative Splicing. RNA-Seq read alignments were further investigated to identify differential splicing between control and ALS samples. Using MATS 24 , a software tools implementing Bayesian multivariate analysis of transcript splicing, we discovered 794 aberrant splicing events (Table 6). Of these, 98% involved exon skipping (including mutually exclusive exons) and 2% altered donor or acceptor splice sites and intron retention ( Table 6).
IPA analysis of gene targets for differential splicing revealed perturbations in genes with important neuronal functions (Supplementary Table 9). In particular, we found genes involved in the actin cytoskeleton signaling, which plays central roles in dynamic processes such as cell motility or axon guidance. GNRH signaling, which coordinates the levels of hormones in the hypothalamic-pituitary-gonadal axis and promotes the activation of G-Proteins, PLC, PKC or the mobilization and influx of calcium was also enriched in genes showing differential splicing as well as genes involved in triacylglycerol biosynthesis and mitochondrial functions (Supplementary Table 9).
Twenty-one of the genes exhibiting differential exon skipping were previously identified in an exon array analysis of motorneurons in sporadic ALS samples 25 .
Gene Ontology terms highly enriched in differentially spliced genes included: axon guidance, neuron projection guidance, anchoring junction and actin binding (Supplementary Table 10). miRNA-Seq analysis. An emerging hallmark of neurodegenerative diseases is the defective RNA metabolism. It involves aberrant alternative splicing as demonstrated also by our RNA-Seq results and deregulated gene expression due to malfunction of miRNA activity 26 . Several previous investigations have documented a crucial role of miRNA gene regulation in the development and function of neurons. Consequently, alterations in the function of miRNA may contribute to neurodegenerative processes. Indeed, specific miRNA alterations have been reported for ALS microglia 27 and affecting the stability of the low molecular weight neurofilament mRNA 28 .
We performed massive sequencing of the low molecular weight RNA fraction extracted from the same control and ALS donor tissues using in the RNA-seq analysis. Between 2.8 and 4.2 million reads were generated per sample and adaptors were removed using a custom script (adapters were detected in between 94.7-95.7% of reads for each sample) (Supplementary Figure 1). Reads were mapped to the human genome and miRNAs (reads mapping perfectly to annotated miRNAs and miRNA* were counted and assigned to non-redundant annotated miRNA sequences).
The vast majority of reads in the peak from 30-34 nt in length corresponded to the 5′ portions of annotated tRNA genes (consistently from nt1 to the anticodon loop) as observed previously 29 . The peak at 18 nt was comprised, predominantly (>95%), of tRNA derived small RNAs, almost exclusively consisting of few fragments  representing positions 1-18 of tRNA-gly and tRNA-ala molecules. None of the aforementioned reads demonstrated statistically significant differential expression between ALS patients and healthy controls. The peak at 21 and 22 nucleotides was strongly enriched from known miRNAs that were subjected to differential expression analysis. After correction for multiple hypothesis testing, we found 18 dysregulated miRNAs (Padj ≤ 0.05) ( Table 7). Of these, 3 were up-regulated and 15 down-regulated in ALS samples.
MiR-155 has previously been shown to be overexpressed in ALS microglial cells 30 and inhibition of its expression has been shown to increase survival time in a murine ALS model 31 . MiR-485 has been implicated in synaptic formation and maintenance and has been reported to be dysregulated in Alzheimer's disease and Huntington's disease 32 , miR-124 is known as a neurodevelopmental regulator 33 , miR-219 required for neural precursor differentiation in zebrafish 34 and promotes myelination in rats 35 , miR-127 regulates cell proliferation and senescence   by targeting BCL6 36 while miR-136, miR-873 and miR-410 have all been studied in the context of glioma growth, invasiveness and apoptosis [37][38][39] . Uncertainties regarding both the reliability of target predictions and the extent of degradation of individual targeted mRNAs in vivo, as well as the mixed nature of the tissue samples complicate direct reconciliation of miRNA and putative target expression profiles. The Diana-miRPath tool 40 recovered a striking set of 6 pathways (p ≤ 10 −16 ) enriched in genes with predicted targets for differentially expressed miRNAs: the PI3K-Akt signaling pathway, Axon guidance, Neurotrophin signaling pathway, Focal adhesion, TGF-beta signaling pathway and Insulin signaling pathway (Supplementary Table 11).

RNA editing analysis.
To further characterize the transcriptome of ALS spinal cord samples, we investigated RNA editing limiting our search to A-to-I changes, as these account for over 99% of editing events in recent large-scale human investigations 41,42 . In particular, for each RNA-Seq sample we explored known A-to-I events stored in the REDIportal database 43 , selecting only genomic positions supported by at least 30 reads and showing a RNA editing level higher than 0.1. Comparing positions in common between healthy and ALS donors, we found a strong correlation of 0,88 (Pval < 10e −14 by Spearman correlation test) in RNA editing levels. We also performed the non-metric multidimensional scaling (NMD) analysis on detected A-to-I changes in order to verify the grouping of input samples by RNA editing. Interestingly, we did not find specific ALS and control clusters as expected in case of global RNA editing alteration (Supplementary Figure 2), suggesting that it should not be impaired in ALS spinal cord tissues.
We observed reduction of RNA editing levels at the recoding GRIA2 Q/R site, known to be edited at nearly 100% in human brain 44 and under-edited in motor neurons of ALS donors and, thus, implicated in cell death by excitotoxicity 45 . However, reduced RNA editing levels at the recoding GRIA2 Q/R site were also observed in control samples and the difference between ALS and control groups was not significant (Mann-Whitney P = 0.43). Only one ALS sample showed strong A-to-I reduction level of 0.19 at GRIA2 Q/R site. Additionally, we compared our RNA editing levels of GRIA2 Q/R site with corresponding values from spinal cord and other brain locations of an independent group of healthy individuals from GTEx project and stored in our REDIportal database. Interestingly, GRIA2 Q/R site appeared highly edited in almost all brain locations (average editing levels from 0.82 to 0.99), while editing levels of spinal cord samples (n = 15, average editing levels of 0.66) were quite similar to our values for both ALS (average editing levels of 0.62) and control (average editing levels of 0.66) groups, meaning that RNA editing levels at GRIA2 Q/R site are not always edited at nearly 100% in spinal cord (Fig. 7A and Supplementary Table 12). Editing frequencies of GRIA2 Q/R site in spinal cord samples were always significantly lower than corresponding values from all tested brain locations (Mann-Whitney P < 0.05) ( Table 8).
Further, we explored RNA editing in miRNAs employing massive sequencing of the low molecular weight RNA. We found 12 A-to-I changes in 12 miRNAs (Supplementary Table 13). Of these, 9 were known as edited in the literature, while 3 appeared as novel RNA editing events. All detected changes were observed in both control and diseased donors, even though RNA editing levels were not significantly different (Supplementary Table 13).
Finally, we explored expression values of ADAR genes but no significant differences were detected between healthy and ALS tissues (Fig. 7B).  Immunoblotting and Immunofluorescence of synaptic SNARE proteins in ventral horns of ALS donors. We corroborated our RNA-Seq and qRT-PCR results performing immunoblotting of two synaptic SNARE proteins, SNAP25 and STX1B, since they might provide new insights into neurotransmission and calcium homeostasis alterations in ALS. Western blot analysis was performed on ventral horn of six patients with sporadic ALS and six controls to quantify the protein expression levels of SNAP25 and STX1B. Results obtained showed that SNAP25 and STX1B were significantly downregulated at protein level in ALS patients, in line with the mRNA data ( Fig. 8 and Supplementary Figure 3). Since a functional redundancy between STX1B and its homologous STX1A has been reported 46 , STX1A protein expression levels were also analyzed in parallel. Results revealed an opposite trend for STX1 A being significantly upregulated in ALS patients ( Fig. 8 and Supplementary  Figure 3), as observed in FUS-silenced motor neurons 47 . Such upregulation appeared more prominent in five out of six ALS patients. Additionally, we investigated in more detail the SNAP25 and STX1A deregulation by immunofluorescence followed by quantification analysis (Figs 9 and 10, Supplementary Figure 4) using confocal microscopy and β3 tubulin as motor neuron marker. The intensity of SNAP25 and STX1A signal was normalized on the size of motor neuron cell body that we found to be reduced by almost 50% in ALS patients (Fig. 9B). In agreement with the SNARE-complex deregulation detected by Western blot (Fig. 8), we found that SNAP25 signal was almost completely absent in ALS patients (Fig. 9), whereas STX1A was significantly upregulated (Fig. 10). Together, these results validate the RNA-Seq data on the altered expression pattern of genes in ALS and provide new information on the molecular basis of synapsis and calcium handling abnormalities associated with the motor neuron degeneration typical of this pathology.

Discussion
ALS is a rare and devastating neurodegenerative disease for which no effective therapies are available 1 . Large-scale genetic investigations, at least in the familial form of ALS, have identified several mutations in key genes for motor neuron survival 2 . Nonetheless, the causes of ALS onset and progression are yet elusive. Whole transcriptome    Table 1. The image in A was cropped for clarity from the original in Supplementary Figure 3.
analyses have revealed the disruption of a variety of cellular pathways contributing to the elucidation of molecular dynamics characterizing the pathobiology of ALS 3 . The main obstacle complicating ALS investigations is that motor neurons cannot be sampled during life. As a consequence, animal models and spinal cord tissue from post-mortem donors represent the elective materials for the study of gene expression in ALS. In humans, previous transcriptome studies have been conducted on laser microdissected motor neurons, in order to avoid biases from mixed cell populations. Although motor neurons are the most affected in ALS, other cell types such as astrocytes, microglia and oligodendrocytes may contribute to disease progression. Indeed, co-culture of motor neurons and glial cells has shown an intricate interplay and an active role of non-neuronal cells in neurodegeneration 6 .
In order to characterize whole transcriptome modifications occurring in sALS, we conducted a study in which total RNA from the lumbar spinal cord of ALS donors were sequenced using the RNA-Seq technology.
Our results demonstrate the suitability of deep transcriptome sequencing for capturing molecular changes in the lumbar spinal cord and provide new exciting insights still not revealed by recent whole transcriptome studies based on RNAseq in brain 7 and cervical spinal cord 8 from fALS and sALS samples.
Using a dataset of cell type-specific genes derived from various cell types in the CNS 15 as astrocytes, oligodendrocytes and neurons, our data show a depletion of neuron, motorneuron and oligodendrocyte populations with down regulated genes of neuronal derivation and linked to relevant neurological functions as impulse transmission, synaptic transmission, calcium ion transport or neurotransmitter secretion. Up regulated genes, instead, are generally involved in neuroinflammation, immune responses or cell death and tend to have microglial origin.
RNA-Seq of ALS lumbar spinal cord therefore provides a snapshot of transcriptomic events characterizing the disease, suggesting that neuronal function is compromised through the down regulation of many genes. As a consequence of neurological injury, surrounding cells are activated leading to a cascade of inflammatory events. Although our data do not provide evidence regarding early molecular events causing the death of motor neurons, transcriptomic results would seem to support these speculations. In addition, comparing our results with those by Brohawn et al. 8 in cervical spinal cord of sALS donors using similar technology and methodology, it emerges that neuroinflammation is common in all sections of spinal cord but molecular signatures of neurodegeneration are more consistent in the lumbar rather than cervical spinal region. Our down regulated genes actually include genes already known to be associated to ALS motor neuron death such as INA, HECW1 or SLC1A2.
As a remarkable result, RNA-Seq data support the down-regulation in ALS donors of SNAP25 and STX1B, neuronal t-SNAREs involved in vesicle trafficking and calcium dynamics. Since a redundant role between STX1B and its homologous STX1A has been reported 48 and being STX1A the isoform interacting with SNAP25 and Voltage Gated Calcium Channels (VGCCs) to control calcium signaling, we have explored the t-SNARE pattern as a whole. In particular, we show that SNAP25 and STX1B are strongly downregulated also at protein level. STX1A, instead, was upregulated in five out of six ALS patients. Through immunofluorescence followed by quantification analysis (Figs 9 and 10) using confocal microscopy and β3 tubulin as motor neuron marker, we confirm the deregulation of SNAP25 and STX1A in ALS patients and demonstrate that the down expression of SNAP25 is mainly due to the reduced size of motor neurons.
Ikemoto et al. 49 have investigated SNAP25 and STX1 expression in ALS spinal cord finding no apparent decrease or only a mild reduction for both proteins in ALS patients. We believe that the immunolocalization approach performed in the study by Ikemoto 49 , without quantification analysis, has most likely led to an underestimation of the t-SNARE reduction in ALS patients. Moreover, the antibody used to detect STX1 was not isoform specific and therefore STX1 unaltered expression reported by Ikemoto 49 is consistent with the data here reported where STX1A is upregulated and STX1B downregulated. Additionally, upregulation of STX1A and deregulation of SNAP25 splicing were observed also in FUS-silenced motor neurons 47 .
Experiments performed in animal models and primary cultures have shown that loss of SNAP25 impairs calcium-evoked exocytosis leaving almost unaffected spontaneous vesicle release 50, 51 while deletion of STX1A has no effect on both types of synaptic transmission probably due to compensation by STX1B. Finally, STX1B deletion impairs both spontaneous and evoked fast synaptic vesicle exocytosis 48 .
According to our results on altered t-SNARE protein pattern, SNAP25 and STX1B reduction seem to strictly relate with the motor neuron dysfunction characteristics of ALS ventral horns in which STX1A up regulation may compensate for STX1B loss. Notably, the observed SNAP25 reduction highlights a potentially novel and relevant role of t-SNARE pattern in elevating intracellular calcium concentration and in triggering the glutamate excitotoxicity. Indeed, SNAP25 silencing in glutamatergic neurons has been reported to increase the concentration of intracellular calcium evoked by depolarization 52 . Therefore an initial accident altering the t-SNARE pattern causing SNAP25 reduction should lead to the increased intracellular calcium levels found in ALS motorneurons 12 . In turn, elevated intracellular calcium concentration causes motor neuron injury leading to increased glutamate release and excitotoxicity. In this context, STX1A may compensate for SNAP25 reduction in order to restore calcium homeostasis. Another remarkable finding of our work is the dysregulation of cholesterol synthesis, which has been linked to chronic neurodegenerative disorders 53 . This alteration could contribute to the impaired axon guidance and synaptic transmission observed in ALS, as cholesterol is essential for the organization and physiological fluidity of cellular membranes.
RNA-Seq also confirms the disruption of RNA metabolism through aberrant alternative splicing potentially resulting from neuronal stress.
Additionally, we investigated RNA editing since its deregulation has been linked to several neurodegenerative diseases including ALS 45 . Indeed, it has been shown that reduced levels of recoding A-to-I editing at GRIA2 Q/R site may be involved in motor neuronal death through an increased Ca 2+ permeability of AMPA receptor 45 . Our results, based on known RNA editing events, do not support significant alterations. A closer investigation of RNA editing levels at the recoding GRIA2 Q/R site reveals that A-to-I frequencies of GRIA2 Q/R site are not impaired in spinal cord tissue of ALS donors, rather their values are similar to those observed in physiological conditions without altered expression of ADAR enzymes. However, our study was conducted on entire tissues and, thus, we cannot exclude a compensating effect from non-neural cells obscuring specific motor neuron RNA editing alterations.
Small RNA expression analyses suggest important roles for miRNAs regulating signaling pathways determining cell survival and axon guidance during ALS progression. Numerous lines of evidence implicate TGF-beta signaling in motor neuron diseases. It participates in the interplay between neuronal and glial cells during ALS progression 6 and increased levels of TGF-beta provoke an acute improvement in the motor performance of SOD1 mice 54 . Participants in the TGF-beta signaling pathway predicted to be targets of DE miRNAs include Smad2,4,5 and 7. TDP-43 and phosphorylated Smad2 are co-localized within cytoplasmic inclusions in the anterior horn cells of sporadic ALS patients 55 . Activation of the TGFβ/Smad signaling system is protective against aggregate formation of cytoplasmically mislocalized TDP-43 55 .
The Axon Guidance pathway has recently been proposed to represent a prime effector in ALS pathobiology 56 and is intimately linked to the regulation of actin cytoskeleton and Focal Adhesion circuits.
The phosphatidylinositol 3-kinase (PI3K) /ATK prosurvival pathway has been implicated as a mediator of the protective effect of vascular endothelial growth factor in G93A-SOD1 neurons in vitro and induction of ATK3 is neuroprotective in SOD1_G93A mice 57 . The pathway contributes, in conjunction with mTOR family proteins, to the regulation of the activity of the glutamate transporter GLT1 in astrocytes 58 .
Neurotrophins, including Nerve Growth Factor (NGF), brain-derived neurotrophic factor (BDNF), and Neurotrophins 3 and 4 (NTF3 and NTF4) lie upstream of the Neurotrophin pathway which regulates the axon cytoskeleton (and in turn axon outgrowth/guidance and synapse formation) as well as cell differentiation and survival. NFT3, NTF4, and BDNF, as well as their cellular receptors and numerous downstream genes are all predicted targets of miRNAs dysregulated in ALS. NGF, BDNF, NT-3 are neuroprotective on axotomized extraocular motoneurons in neonatal rats 59 , while 7,8-dihydroxyflavone (7,8-DHF) a small molecule tyrosine kinase receptor B (TrkB) agonist that mimics the effects of BDNF, is neuroprotective in SOD1 mice 60 .
In conclusion, our study remarks the use of massive RNA sequencing technologies to improve deregulated pathways of ALS neurodegeneration and provide new exiting insights as the SNAP25 reduction as possible cause of calcium elevation and glutamate excitotoxicity, suggesting t-SNARE protein expression as novel indicators and potential biomarkers for sporadic ALS.

Methods
Post mortem tissues. Frozen lumbar spinal cord samples from 11 male human donors affected by sporadic ALS and 7 age, sex and ethnicity matched controls were obtained from the NICHD Brain & Tissue Bank for Developmental Disorders (University of Maryland -http://medschool.umaryland.edu/btbank/) and the Human Brain and Spinal Fluid Resource Center (Los Angeles, CA). Additionally, we required as an additional criterion, the presence of clinical symptoms associated to the altered function of lower motor neurons.
Each tissue was longitudinally dissected and ventral horns were collected for downstream analyses. Detailed information regarding samples used in RNA-seq, miRNA-seq, qRT-PCR and western blot analyses are provided in Table 1, also including age, gender, ethnicity and post-mortem interval as well as donor bank IDs to access to further clinical info. Preprocessing and analysis of RNA-Seq reads. RNA-Seq reads in FASTQ format were initially inspected using FASTQC program (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Adaptors and low quality regions (phred cutoff of 20) were trimmed using Trim Galore (http://www.bioinformatics.babraham. ac.uk/projects/trim_galore/), excluding reads with final length less than 50 bases.
Cleaned reads were aligned onto the complete human genome (assembly hg19) by means of GSNAP 61 version 2013-11-27 (using as parameters: -B 5 -d hg19 -t5 -s splicesites -E1000 -N1 -n1 -Q -O-nofails -A sam-force-xs-dir -a paired) providing a list of exon-exon junctions from Ensembl, UCSC and RefSeq databases. Unique and concordant alignments in SAM format were converted in the binary BAM format by SAMtools 62 and basic statistics were calculated using the CollectRnaSeqMetrics tool of Picard package (https://broadinstitute.github.io/picard/ index.html) ( Table 2). The sporadic disease nature of our samples was checked using REDItools 13 , screening known ALS mutations from ALSoD database in aligned RNAseq data 14 . Only positions showing a frequency variation higher than 0.1 and supported by more than 10 reads we taken into account.
Differential expression was performed using two independent tools: CuffDiff2 version 2.1.1 17 (using as main parameters:-library-type fr-firststrand-labels Ctl,Als -u -b hg19.fa refgenes.gtf) and DESeq2 version 1.10.1 18 . Read counts on known human genes were calculated by featureCounts version 1.5.1 63 . Reference human transcriptome was obtained from iGenomes repository (http://cole-trapnell-lab.github.io/cufflinks/ igenome_table/) and annotations for rRNA genes were downloaded from UCSC genome browser selecting the RepeatMask table. DE genes was selected intersecting lists of significant DE genes (corrected P value < 0.05) obtained by CuffDiff2 and DEseq2. We took into account only genes with |log2(FC)| higher than 1.
Changes in cell type composition of spinal cord samples was detected according to the method described in Lin et al. 64 and using cell type specific genes for neurons, motorneurons, astrocytes, oligodendrocytes and microglia collected from the literature 6,15 . The majority of cell type specific genes were obtained from Cahoy et al. 15 taking into account only genes showing >10-fold enrichment in this study. Since Cahoy et al. performed their investigations in mouse, we selected only genes with corresponding human orthologues and showing FPKM values > 0. Motor neuron specific genes, instead, were collected from Phatnani et al. 6 .
Differential alternative splicing was detected with MATS version 3.0.9 24 using default parameters except for the read length (fixed to 90). Events were selected at 0.05 significant level, corrected for multiple testing. Heat-maps and PCA plots were generated using ClustVis 65 .
IPA system was used for pathway analysis. Enrich 66 , instead, was used to generate lists of Gene Ontology terms enriched in the upregulated, downregulated or differentially spliced gene lists, with default setting. To identify sets of genes that are both strongly correlated to expression data and functionally related to their GO annotations, we used the GO-PCA method 23 . In brief, it adopts a two-step approach in which PCA is performed first. Then, each principal component is tested for whether it is driven by functionally related genes (in the form of gene ontology annotations).
Preprocessing and analysis of miRNA-Seq reads. Raw sequence data were processed with a custom Python script to remove adapters from inserts of between 14 and 38 nt in length and subsequently mapped onto the human genome (version hg19) using Bowtie 67 . Counts of reads exactly matching human miRNAs and (where annotated in miRBase) miRNA* sequences were submitted to the DESeq software 68 for statistical analysis of differential expression using default parameters and treating different ALS and healthy control patients as biological replicates. Pathway analysis was conducted with miRPath v2.0 40 with default settings.
RNA editing detection and analysis. RNA editing candidates per each sample were detected using the REDItools package (https://sourceforge.net/projects/reditools/) 13 . In particular, REDItoolKnown.py was used to explore known RNA editing sites stored in REDIPortal database 43 . Non-Metric Multidimensional Scaling analysis of RNA editing sites was performed in R using the vegan package. RNA editing levels at the GRIA2 Q/R site in different brain regions were extracted from REDIportal 43 . Custom python scripts were used to calculate RNA editing correlation between control and ALS samples.
RNA editing detection in miRNAs was performed according to the protocol developed by Alon et al. 69 , while RNA editing levels were calculated by REDItools 13 .
qRT-PCR validation. Gene expression analyses were carried out on custom TaqMan Array Fast Plates (Life Technologies), configured in order to include four potential human endogenous control genes (GAPDH, HPRT1, PGK1 and POLR2A), and sixteen genes selected for validation ( Table 5), each of them in duplicate per sample. Reverse transcription of 500 ng of total RNA was performed using "iScript Reverse Transcription Supermix for RT-qPCR" (Bio-Rad, Hercules, CA, USA). qPCR reactions were set up according to the manufacturer's instruction of "TaqMan Fast Universal PCR Master Mix" kit (Life Technologies), starting from 1 µl of diluted cDNA (1:3) as template and performed using the ABI PRISM 7900HT platform (Applied Biosystem, Life Technologies) with the following amplification conditions: hot start at 95 °C for 10 min; 45 amplification cycles (95° for 15 sec, 60° for 1 min).
Fluorescence raw data were exported from the SDS 2.2.1 software and analyzed using the comparative quantification cycle (Cq) method to determine the relative ratio of transcripts. First, the ΔCt values were obtained by calculating, for each sample, the difference of the Cq of each target compared to the arithmetic mean of the endogenous Cq genes (GAPDH and HPRT1). Then, the relative expression ratio (rER) was expressed as SCIENtIfIC REPORtS | 7: 10046 | DOI:10.1038/s41598-017-10488-7 Log2(2 −ΔCt ); rER for each type of samples were averaged and calibrated respect to normal samples. Analyses were carried out using data related to two independent experiments.
The following secondary antibodies were used for Western blot: HorseRadishPeroxidase (HRP)-conjugated donkey anti-goat, goat anti-rabbit and goat anti-mouse IgG (Santa Cruz Biotechnology, Santa Cruz, CA, USA).
Protein sample preparation and immunoblotting. 1 mm slices of Ctrl and ALS ventral horn samples, were collected, cut into small pieces with a razor blade and immediately solubilized in 8 volumes of RIPA buffer (10 mM Tris-HCl, pH 7.4; 140 mM NaCl; 1% Triton X-100; 1% Na deoxycholate; 0.1% SDS; 1 mM Na3VO4; 1 mM NaF and 1 mM EDTA) added with a cocktail of protease inhibitors (Roche, Milan, Italy). The lysis was performed as previously described 70 . Briefly, on ice for 1 h and the samples were then centrifuged at 22,000xg for 1 h. The protein content of the supernatant was measured with a bicinchoninic acid (BCA) Protein Assay Kit (Rockford, IL, USA). 30 ug of protein samples were separated by 13% Tris-Glycine-SDS-PAGE and transferred to polyvinylidene difluoride membranes (Millipore, Milan, Italy) as described previously 71 . Membranes with blotted proteins were incubated with primary antibodies, washed, and incubated with peroxidase-conjugated secondary antibodies. Reactive proteins were revealed with an enhanced chemiluminescent detection system (ECL Plus; GE Healthcare, Buckinghamshire, UK) and visualized on a Versadoc imaging system (BioRad, Milan, Italy). Densitometry analysis was performed using Scion Image software (Frederick, MD, USA).
All data are reported as mean ± SEM. We used the Student's t test for unpaired data and analysis of variance for multiple statistical comparisons between groups, the significance level being set at P < 0.05. The number of samples is indicated in the figure legend.
Immunofluorescence. Fresh sections of control and ALS spinal cord samples were rehydrated with PBS and fixed for 10 min with 4% paraformaldehyde in PBS (HT5014, Sigma, Milan Italy). Sections were washed in PBS and incubated with 1:200 diluted anti-SNAP25 and 1:1000 diluted anti-β3 tubulin or 1:200 diluted anti-Syntaxin-1A and 1:1000 diluted anti-β3 tubulin in 0.3% Triton X-100 in PBS overnight at 4 °C. After washings with PBS, sections were saturated with 0.1% Gelatin (Sigma, Milan, Italy) in PBS, and incubated with 1:1000 diluted AlexaFluor488 and AlexaFluor594 conjugated secondary antibodies for 1 hour at room temperature. After washings, sections were washed in PBS and mounted with Mowiol (Sima, Milan, Italy) and acquired using Leica TCS SP2 confocal microscope (Leica Microsystems, Heidelberg GmbH) for the quantitative analysis described in the next paragraph. Confocal 3D reconstruction was obtained using Leica TCS SP8 STED 3X and Leica HC PL APO 100x/1.40 Oil white objective with Type F Immersion liquid with 1.5 refractive index (Pisani et al., 2016). Excitation of AlexaFluor488 and AlexaFluor594 dyes was performed with a continuous 488 nm and 594 nm wavelength diode laser with a maxim light output in focal plane of 10 mW (NKT Photonics sign supercontinuum laser). To obtain the confocal 3D projection, 30 planes were acquired in a Z-stack of 6-12 μm, and the series processed by Leica LASX software (Leica Microsystems, Mannheim, Germany).
Quantitative analysis of motor neuron size and SNAP25 and STX1A expression by immunofluorescence in ventral horn sections. Single confocal planes of randomly chosen areas of ventral horn sections from three control donors and three ALS patients, stained as described above, were acquired by Leica TCS SP2 confocal microscope (Leica Microsystems, Heidelberg GmbH). Leica Confocal Software Version 2.61 was used to measure motor neuron size and to quantify SNAP25 and STX1A expression levels. Motor neuron cell bodies were drawn based on the β3 tubulin staining, and the surface quantified and expressed in µm 2 . The average signal value relative to SNAP25 and STX1A for each motor neuron was also calculated and normalized on the motor neuron surface. Results were expressed as means ± standard error. Statistical significance was evaluated using unpaired Student's t test. Difference with p < 0.05 were considered statistically significant. Data availability. All sequencing data produced in the present work have been submitted to dbGaP database under the accession phs000747.