Elucidating the Contribution of Skeletal Muscle Ion Channels to Amyotrophic Lateral Sclerosis in search of new therapeutic options

The discovery of pathogenetic mechanisms is essential to identify new therapeutic approaches in Amyotrophic Lateral Sclerosis (ALS). Here we investigated the role of the most important ion channels in skeletal muscle of an ALS animal model (MLC/SOD1G93A) carrying a mutated SOD1 exclusively in this tissue, avoiding motor-neuron involvement. Ion channels are fundamental proteins for muscle function, and also to sustain neuromuscular junction and nerve integrity. By a multivariate statistical analysis, using machine learning algorithms, we identified the discriminant genes in MLC/SOD1G93A mice. Surprisingly, the expression of ClC-1 chloride channel, present only in skeletal muscle, was reduced. Also, the expression of Protein Kinase-C, known to control ClC-1 activity, was increased, causing its inhibition. The functional characterization confirmed the reduction of ClC-1 activity, leading to hyperexcitability and impaired relaxation. The increased expression of ion channel coupled AMPA-receptor may contribute to sustained depolarization and functional impairment. Also, the decreased expression of irisin, a muscle-secreted peptide protecting brain function, may disturb muscle-nerve connection. Interestingly, the in-vitro application of chelerythrine or acetazolamide, restored ClC-1 activity and sarcolemma hyperexcitability in these mice. These findings show that ion channel function impairment in skeletal muscle may lead to motor-neuron increased vulnerability, and opens the possibility to investigate on new compounds as promising therapy.


Results
Gene expression profile in skeletal muscle of SOD1 G93A and MLC/SOD1 G93A mice. The quantitative Real-Time PCR analysis have been used to detect skeletal muscle alterations due to SOD1 mutation in both SOD1 G93A and MLC/SOD1 G93A mice with respect to age-matched and strain-matched wild-type (WT). We evaluated the expression of a series of genes involved in skeletal muscle structure and function, encoding for ion channels, their regulatory proteins and subunits, structural proteins, markers of denervation ( Supplementary  Fig. S1) and phenotype shift ( Supplementary Fig. S2), hormones and neurotrophic factors. The involvement of these genes in ALS etiology was evaluated by an accurate multivariate statistical analysis PCA-LDA (Principal Component Analysis and Linear Discriminant Analysis), using both unsupervised and supervised machine learning algorithms, which allowed to identify new potential therapeutic targets.
Multivariate Statistical Analysis of gene expression data. Principal Component Analysis (PCA) and Linear Discriminant Analysis (LDA) The gene expression dataset was composed of 27 columns representing the expression level of genes normalized to β-Actin (housekeeping gene) and 30 rows representing the analyzed TA muscles. This dataset has been used to perform a PCA calculating 27 eigenpairs of eingenvector-eigenvalue to project in the newly modelled subspace of characteristics. This analysis shows that the first 3 Principal Components (PCs) accounted for about the 2/3 of the total of the variance ( Supplementary Fig. S3, Tables S1 and S2). Muscles of SOD1 G93A group are well separated from the others ( Fig. 2A-C). The selected genes were able to separate well also the less affected MLC/SOD1 G93A group from their WT. These observations were evident after analysis of the scatterplots in Fig. 2D. In order to confirm it mathematically, we used the PCA-LDA approach. Thus, training an LDA algorithm, we obtained the contribution of each PC to separate the three pairs of groups: (1) SOD1 G93A and their age-and strain-matched controls; (2) MLC/SOD1 G93A and their strain-matched controls; (3) MLC/SOD1 G93A and SOD1 G93A . The three LDA (Fig. 2E) were able to find directions (Linear Discriminant, LD) that well segregate samples. The LD coefficients indicated the importance of each PC to separate groups (Supplementary  Table S3). Through linear combination of LD coefficients and the PC weight of each gene transcripts (PCA-LDA www.nature.com/scientificreports www.nature.com/scientificreports/ pipeline), we computed three lists of discriminant genes (Fig. 3, Supplementary Table 4) showing a high weight in the distinction between groups. Farther from the zero the weight is, the greater its discriminating ability. For the separation of MLC/SOD1 G93A and WT, the genes with the higher discriminant power were AMPAR2, PKC-theta, Irisin. For MLC/SOD1 G93A and SOD1 G93A pairs, the highest discriminant genes were Myosin Heavy Chain-2b (MyHC-2B), Acetylcholine Receptor Subunit alpha-1 (AChRa1), SERCA2, ClC1, Myosin Heavy Chain-2a (MyHC-2A), Nav1.4. For separation between SOD1 G93A and WT, the following genes were discriminant: AChRa1, SERCA2, MyHC-2B, MyHC-2A, Notch1, Sur1, Myosin Heavy Chain-1 (MyHC-1), ClC-1.

Hierarchical agglomerative clustering.
To confirm the pattern modification of transcript levels, we used a hierarchical agglomerative clustering algorithm (Fig. 4). This analysis showed, in an unsupervised manner, that both genes and muscles segregated into two big different clusters. Also in this case, the clusters segregated muscles belonging to ubiquitously SOD1 G93A into (MLC/SOD1 G93A ) and WT, highlighting that the major part of the variance is due to variation in gene expression level occurring in SOD1 G93A mice.
ClC-1 protein expression measured by western blot and immunofluorescence assay in SOD1 G93A and MLC/SOD1 G93A mice. To better explore on the ClC-1 channel modification in skeletal muscle of ALS mice we analyzed the expression of ClC-1 protein by Western Blot (WB) analysis in TA muscles. We found a significant reduction of ClC-1 expression with respect to WT in both models (Fig. 5A), suggesting that a significant less amount of channel protein contributes to the gCl. Since the mRNA level of ClC-1 gene was not significantly changed in the muscle of MLC/SOD1 G93A mice, it seems that post-transcriptional events are involved, i.e. an alteration of translation or a rapid degradation of the protein. We also evaluated, by an immunofluorescence assay, the presence of the ClC-1 in cryosections of TA muscle excised from MLC/SOD1 G93A mice and of gastrocnemius muscle from SOD1 G93A mice (Fig. 5B) both with high expression level of the transgene 2 . The presence of ClC1 protein in muscle sections was shown by using specific antibodies and no modifications were observed in both situations (Fig. 5B). We measured the resting component conductances (gCl and gK: Potassium Conductance) in EDL muscle of both ALS mouse models. In detail, SOD1 G93A muscle fibers showed a significant reduction of gCl with respect to their strain-and age-matched controls by −27.8 ± 4.2% and −33.9 ± 3.8%, at 90 and 130 days of age, respectively (Fig. 5C,D). Thus, a reduction of this parameter was already present at clinical onset, and was more severe at the exacerbation of the signs of pathology. A slight increase of gK (by +32.8 ± 11.9%) was found in 90 days-old SOD1 G93A animals with respect to their controls, while the difference was significant (by +67.5 ± 27.9%) www.nature.com/scientificreports www.nature.com/scientificreports/ in 130 days-old SOD1 G93A mice (Fig. 5F). In MLC/SOD1 G93A mice, we found a significant reduction of gCl by −22.3 ± 3.2% with respect to the WT (Fig. 5E). However, no significant differences of gK were observed (Fig. 5G). No modification was found in gCl and gK measured in slow-twitch SOLEUS (SOL) muscles of SOD1 G93A and WT at 130 days of age ( Supplementary Fig. S4). To evaluate the involvement of the PKC in the reduction of gCl we tested the effect of chelerythrine, a selective inhibitor of PKC 23 . We found a significant restoration of gCl in EDL fibers of SOD1 G93A mice at both ages and of MLC/SOD1 G93A mice after in vitro application of chelerythrine  . Hierarchical agglomerative clustering of gene expression data. Heatmap of unsupervised machine learning algorithm to perform clusterization of muscles based on gene expression levels. TA muscles and mRNA levels segregated into two principal clusters also in this case, as showed by hierarchical dendrograms of rows and columns. Each square represents the standardized gene expression level normalized with β-actin housekeeping gene in TA muscle of target genes. The columns report target genes, rows report the analysed muscles. Genes and muscles are sorted by hierarchical agglomerative clustering algorithm. Plotted data consists of the entire dataset of animals (SOD1 G93A , MLC/SOD1 G93A and their strain-matched WT as control). (C,D) The macroscopic Chloride Conductance (gCl) measured in EDL muscles fibers of SOD1 G93A animals at the two selected ages and effects of the in-vitro application of chelerythrine (1 µM). Statistical analysis was performed using one-way ANOVA followed by Bonferroni post-hoc t-test (F = 15.8, df = 2, P < 0.0001, for gCl in SOD1 G93A at 90-days) (F = 16.6, df = 2, P < 0.0001 for gCl in SOD1 G93A at 130-days). *Significantly different vs. age-matched WT (at least P < 0.05). (E) Macroscopic gCl measured in EDL muscles fibers of MLC/ www.nature.com/scientificreports www.nature.com/scientificreports/ ( Fig. 5C-E). In addition, we tested the effects of acetazolamide (ACTZ), an inhibitor of carbonic anhydrase, which has been previously shown to increase heterologously-expressed ClC-1 channel activity 11,14 . Importantly, the in vitro application of ACTZ led to an increase of gCl in EDL muscle fibers of MLC/SOD1 G93A mice toward the WT value ( Fig. 5E) suggesting potential beneficial effects in ALS models characterized by gCl reduction and hyperexcitability.
Excitability parameters in EDL muscle of SOD1 G93A and MLC/SOD1 G93A mice. The excitability parameters were recorded in EDL muscle of SOD1 G93A mice and in muscle specific MLC/SOD1 G93A transgenic mice. Representative traces of the active parameters measured in EDL muscle of SOD1 G93A mice and relative WT, are shown in Fig. 6(A-D). In accord with the reduction of gCl the excitability parameters were modified ( Fig. 6E-H). In particular, the Latency of action potential (Lat, the delay from the beginning of the current pulse to the onset of an action potential at threshold, in milliseconds) was increased in SOD1 G93A mice at 90 and at 130 days of age (+117.7 ± 29.9% and +87.7 ± 20.8%, respectively) with respect to their age-matched WT. The threshold current for the onset of the first action potential (Ith, in nanoamperes) was decreased at both age-points (−64.8 ± 5.6% and −44.4 ± 9.6%, respectively) and the maximum number of elicited spikes was increased (+55.4 ± 15% and +77.5 ± 18%, respectively). Also, we found a decrease of the amplitude (AP, in millivolts) of action potential (by −22.1 ± 6.3%) only at 130 days of age. The excitability parameters of sarcolemma have been measured also in EDL muscle fibers of MLC/SOD1 G93A mice and their strain-matched controls (Fig. 6I). In general, active parameters of MLC/SOD1 G93A mice were less affected with respect to SOD1 G93A ones. Indeed, only the maximum number of elicited spikes was significantly increased by +72.3 ± 11.7% with respect to WT. Importantly, the application of 50 µM ACTZ significantly restored the maximum number of spikes toward control value.

ATP-sensitive potassium channels (KATP) activity in skeletal muscle of SOD1 G93A mice.
Patch-clamp experiments were performed on fast-twitch flexor digitorum brevis (FDB) muscle fibers to test whether the ubiquitous expression of SOD1 G93A can modify the activity of KATP channels, which are known to be altered during various physiopathological conditions 24 . The current amplitude of the KATP channel was increased in the FDB muscle of SOD1 G93A mice with respect to WT at 130 days (Fig. 7A). Exposure of macropatches to intracellular ATP (100µM-5mM) greatly reduced the current amplitude in control conditions, confirming that these currents flowed through KATP channels (Fig. 7B). However, KATP channels of SOD1 G93A mice were less sensitive to 100 µM ATP.
Calcium Homeostasis and response to caffeine in skeletal muscle of SOD1 G93A and of MLC/ SOD1 G93A mice. Using the FURA-2 cytofluorimetric technique, we measured calcium homeostasis in EDL muscle fibers of SOD1 G93A and MLC/SOD1 G93A mice in comparison to strain-matched WT (Fig. 8,B). Muscle fibers of SOD1 G93A mice showed a significant increase in resting cytosolic calcium (restCa) concentration with respect to WT. In-vitro application of 40 mM caffeine, a modulator of RyR, on EDL fibers of SOD1 G93A mice, induced a slight increase of restCa suggesting higher sensitivity to the compound (Fig. 8C,D). The restCa level was not significantly modified in skeletal muscle of MLC/SOD1 G93A mice.

Discussion
In this study we tested for the first time the activity and expression of various skeletal muscle ion channels, which are pivotal for the function of muscle and of neighboring tissues. The analytical strategy performed allowed the identification of skeletal muscle expression signature, based both on differential gene expression and on gene correlation networks. The results show a clear-cut segregation of muscles in separated groups: SOD1 G93A and MLC/SOD1 G93A versus their strain-matched WT. The PCA analysis of gene expression data showed no separation between 90 days-old and 130 days-old animals suggesting that skeletal muscle is severely involved in ALS pathogenesis already before the appearance of clear symptoms. In ALS patients, muscle denervation and abnormal reinnervation are pathological hallmarks of the disease. Similarly, an up-regulation of denervation markers (Nav1.5, AChR) 8,25,26 was found in SOD1 G93A muscles at both time points, together with an increase of myogenin 27-29 that indicates muscle reinnervation process 27 . In addition, inhibition of Mstn, a negative regulator of muscle mass, suggests an attempt of the muscle to counteract atrophy and promote regeneration. Interestingly, Mstn expression was decreased also in MLC/SOD1 G93A mice, the animal model expressing the SOD1 mutant gene selectively in skeletal muscle 4 , again suggesting the early involvement of skeletal muscle in the pathogenesis of ALS. A fast-to-slow phenotype transition was observed in SOD1 G93A mice, likely as an attempt to convert muscle fibers into a more resistant slow type, less affected by the pathology. In skeletal muscles of these mice, the expression of ClC-1 channel and the related resting gCl were significantly decreased with respect to WT, at 90 and 130 days of life. This is in accord with the fast-to-slow phenotype transition, since ClC-1 expression and resting gCl in slow-twitch muscles are lower than in the fast ones 30 . The ClC-1 protein expression and resting gCl were SOD1 G93A animals and effect of in-vitro application of Chelerythrine (1 µM) or Acetazolamide (ACTZ, 50 µM). Statistical analysis was performed using one-way ANOVA followed by Bonferroni post-hoc t-test (F = 12.5, df = 3, P < 0.0001, for gCl in MLC/SOD1 G93A ). *Significantly different vs. age-matched WT (at least P < 0.05). °Significantly different vs. MLC/SOD1 G93A (at least P < 0.05). (F) The macroscopic Potassium Conductance (gK) measured in EDL muscles fibers of SOD1 G93A animals at the two selected ages (F = 5.18, df = 3, P < 0.005, for gK in SOD1 G93A ). *Significantly different vs. age-matched WT (at least P < 0.05). (G) The macroscopic gK measured in EDL muscle of MLC/SOD1 G93A mice (no significant differences were found). (2019) 9:3185 | https://doi.org/10.1038/s41598-019-39676-3 www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/ significantly reduced also in muscles of MLC/SOD1 G93A mice, likely through a mechanism involving Reactive Oxygen Species (ROS) production. An in-depth analysis revealed that the reduction of gCl in both models also involved PKC-theta up-regulation, which phosphorylates and inhibits the ClC-1 channel 13,31 . According to the reduction of gCl, sarcolemma excitability was increased, with a significant increase of the maximum number of spikes. The situation was more complex in SOD1 G93A mice, since other channels were involved. In these mice a decreased amplitude of the action potential observed in SOD1 G93A mice can be in part explained by the down-regulation of Nav1.4 sodium channel transcript. The increase in resting gK is likely sustained by the increased KATP current in SOD1 G93A mice, although we cannot exclude the contribution of other K channels. The KATP current is myoprotective and may be activated to dampen excitability. However, the mRNA expression of the pore subunit Kir6.2 and the auxiliary subunit SUR2A was decreased, but it is possible that the increased expression of SUR1, more resistant to oxidative stress, may account for the increased KATP current. Moreover, the modification of the channel and/or the phenotype shift 32 may explain the ATP resistance in these mice. All these modifications were not observed in the MLC/SOD1 G93A mice. Likewise, calcium channels were strongly involved in the modification of skeletal muscle function in SOD1 G93A mice. The restCa was significantly increased in the EDL muscles according to the fast-to-slow phenotype transition. Despite a reduction of RyR mRNA expression, a strong decrease of SERCA1, responsible for SR calcium reuptake, may contribute to the increase of restCa. Also a depletion of ATP may reduce SERCA1 activity. Interestingly, the increase of SERCA2 expression may represent an attempt to compensate for the loss of SERCA1 [33][34][35] . In MLC/SOD1 G93A model, the not significant increase of restCa and the unchanged expression of RyR and SERCA1, both suggest a minor involvement in this strain. Other studies 36 have shown an impairment of EC coupling in FDB fast-twitch fibers of both models. All the modifications observed in skeletal muscle of SOD1 G93A animals seems to be associated to denervation and phenotype transition. In contrast, no sign of morphological denervation was found in MLC/SOD1 G93A mice, at the examined age, where motor neurons are not affected. However muscle restricted expression of SOD1 G93A displays alteration in NMJ 4 . MLC/SOD1 G93A muscles were characterized by a minor number of transcriptional modifications, and the identified most discriminant genes involved in the onset of damage were PKC-theta, that control muscle ClC-1 activity, AMPAR2, the ionotropic glutamate receptor involved in synaptic plasticity, and irisin, a recently-discovered hormone secreted by skeletal muscle and able to communicate in a paracrine and endocrine manner with nervous system 37 . The functional significance of an increase of AMPAR2 may be the facilitation of glutamate-induced prolonged depolarization and modulation of messengers (e.g. NO) able to activate retrograde potentiation in pre-synaptic compartments 38,39 . The consequent muscle hyperexcitability may rapidly exacerbate the damage at neuronal tissue 9 , already compromised by the increase of firing rates, ionic imbalance and sustained depolarization that decreases ATP level 21,35 . The results indicate that muscle hyperexcitability is due to ClC-1 expression reduction and increased activity of PKC-theta with consequent gCl reduction in MLC/ SOD1 G93A mice. The increased activity of PKC was confirmed by the in-vitro effects of chelerythrine, a well known inhibitor of PKC, which was able to increase the resting gCl 23 . Also, AMPAR2 activation can be modulated by PKC-induced increased phosphorylation 40 . Our hypothesis is that selective accumulation of mutant SOD1 in muscle and consequent oxidative stress may affect PKC, ion channel function, and neuromuscular communication 41 . In this situation, a "feedback-loop" could be triggered, which may explain why motor neuron death begins with an initial distal hyperexcitability [42][43][44] . Accordingly, muscle-specific expression of local Igf-1 (mIgf-1) was shown to stabilize NMJ, reduce inflammation in the spinal cord, and enhance motor neuronal survival in SOD1 G93A mice, thereby delaying the onset and progression of the disease 5,8 . We already demonstrated that Igf-1 signaling is important to maintain a normal phosphatase activity that counterbalances the inhibitory activity of PKC on ClC-1 45,46 . Thus, it is likely that the decrease of phosphatase activity together with PKC activation strongly contribute to gCl reduction and hyperexcitability in ALS. Since gCl is strongly reduced during aging in skeletal muscle of rodents 12 , it is possible that more severe alterations in MLC/SOD1 G93A mice may occur at advanced age, as already observed 5 . Accordingly, preliminary data showed a more severe reduction of gCl in 10-months-old MLC/SOD1 G93A mice (Supplementary Fig. S6). All these findings identify ACTZ as a potential drug for ALS therapy. As already demonstrated 11,14,47,48 , this carbonic anhydrase inhibitor, was able to increase ClC-1 channel activity likely through modification of internal pH. The increased expression of carbonic anhydrase in motor neurons of ALS patients also supports ACTZ therapeutic proposal 49 . Notwithstanding the encouraging results in the restoration of gCl and attenuation of muscle hyperexcitability in MLC/SOD1 G93A mice, additional preclinical studies are required to better assess the long-term effects of ACTZ. Irisin is another discriminant gene, found to be down-regulated in MLC/SOD1 G93A mice. Irisin is a myokine secreted by both skeletal muscle and brain. Physical activity modifies plasma concentration of this hormone known to stimulate the synthesis of brain-derived neurotrophic factor (BDNF) in the nervous system 50 . Thus, irisin associates skeletal muscle activity to that of nervous system in a centripetal manner, being able to cross the blood-brain barrier 51 . Irisin gene expression was modified not only in the MLC/SOD1 G93A mice, but also in SOD1 G93A mice, suggesting that its modification is an early and long-lasting event in the pathogenesis of ALS. Based on these results, irisin may represent a critical link between muscle and CNS in ALS and a likely pharmacological target. Thus, a combined therapy with irisin may be also considered to restore muscle-nerve communication. ACTZ, which increases ClC-1 activity, and riluzole, which mean ± S.E.M. from 10-20 fibers. The excitability parameters were recorded in the absence and in the presence of 50 µM ACTZ. To note, because of the Log scale, the measure units are omitted and are: for Action Potential amplitude (AP) mV; for the threshold current amplitude (Ith) nA; for the latency of AP (Lat) msec; and for the maximum number of elicitable AP (N spikes). Statistical analysis was performed using one-way ANOVA followed by Bonferroni post-hoc t-test (F = 26.86, df = 2, P < 0.0001, for N spikes). *Significantly different vs. WT (P < 0.05) °Significantly different vs. MLC/SOD1 G93A (P < 0.05). (2019) 9:3185 | https://doi.org/10.1038/s41598-019-39676-3 www.nature.com/scientificreports www.nature.com/scientificreports/ inhibit Na + channels activity 14 , might act synergically to restore muscle excitability. Our results strengthen the evidence for the role of skeletal muscle in ALS pathogenesis and pave the way for the development of new therapeutic options to hamper the clinical effects of the disease.

Material and Methods
Animal care. In this study two different animal models of ALS were evaluated. Male SOD1 G93A transgenic mice ubiquitously expressing the human mutant SOD1 G93A allele containing the Gly93RAla (G93A) substitution were used at two time points (90 and 130 days-old) based on the criteria for disease progression 52 as clinical onset (mice selected showed no overt sign of pathology) and early symptomatic (mice were selected at abnormal gait and before paralysis). Male MLC/SOD1 G93A mice (140 days-old) overexpressing the mutant SOD1G93A transgene under the transcriptional control of the Myosin Light Chain (MLC) muscle specific promoter, so that its expression is selectively restricted in skeletal muscle 4 . Transgenic animals were identified by PCR amplification of DNA extracted from tail tissue 44 . Age-matched C57BL6J and FVB/NJ mice (Jackson Laboratories, Bar Harbor, ME, USA) were used as WT control strain, respectively. The animals were housed in a temperature-controlled (22 °C) room with a 12:12 h light-dark cycle. Before muscle dissection animals were deeply anesthetized with ketamine (100 mg/kg ip)/xylazine ( Isolation of total RNA, reverse transcription and real-time quantitative polymerase chain reaction (pCR) analysis. Immediately after surgery, TA muscles were snap frozen in liquid nitrogen and stored at −80 °C until isolation of total RNA, reverse transcription and real-time PCR analysis, as detailed in Camerino et al. 16 . Genes were analysed by the use of TaqMan Hydrolysis primer and probe gene expression assays that are produced by Life-Technologies with the following assay IDs: PKC theta (encoded by Prkcq gene) assay ID: Mm00436796_ m1; PKC alpha (encoded by Prkca gene) assay ID: Mm00440858_m1; Kir 6.2 (encoded by Kcnj11 gene) assay ID: Mm00440050_g1; Sur1 (encoded by Abcc8 gene) assay ID: mm008003450_m1; Ryr1 (encoded by Ryr1 gene) Figure 7. KATP channel current density from WT and SOD1 G93A mice measured by patch clamp technique. Each bar is the mean ± SEM from the number of patches indicated in brackets. (A) significant difference among groups was found by one-way ANOVA analysis (F = 7.2; P < 0.0002). Bonferroni post hoc correction for individual differences between groups is as follows: significantly different * vs WT 130 days (P < 0.002); # vs SOD1 G93A 90 days (P < 0.002). (B) Exposure of macropatches to intracellular ATP (100 µM-5 mM) induce a KATP current inhibition. The response of KATP channels to 100 µM ATP is changed in 130 days-old SOD1 G93A mice.  53 . The mRNA expression of the genes was normalized to the best housekeeping gene: Actb selected from B2m and Actb by Normfinder software 16 . For genes that are poorly expressed, such as Ampar2, Nmdar1, Notch1, Ngf, Mstn pre-amplification by TaqMan PreAmp Master Mix (Life Technologies C.N. 4391128) was made before the real-time experiments.
Western Blot analysis. ClC1 protein was isolated from TA muscle of MLC/SOD1 G93A and processed as described in Camerino et al. 16 . Primary rabbit anti-ClC1 antibody (MyBiosource) was diluted 1:200 with TBS containing 5% non-fat dry milk overnight. Membranes were incubated for 1 h with secondary antibody labeled with peroxidase (1:5000 anti rabbit IgG, Sigma-Aldrich), developed with a chemiluminescent substrate (Clarity Western ECL Substrate, Bio-Rad) and visualized on a Chemidoc imaging system (Bio-Rad). Densitometric analysis of each experimental band was performed using Image Lab software (Bio-Rad). The software allows the chemiluminescence detection of each experimental protein band to obtain the absolute signal intensity. The density volume was automatically adjusted by subtracting the local background. For each sample, the relative intensity was calculated by normalizing the intensity of β-actin protein band (diluted 1:300, rabbit Anti Actin, Santa Cruz Biotechnology) as reference standard.
Immunofluorescence. TA muscles, covered with tissue-tek O.C.T. (Bio-Optica), were frozen in isopentane cooled in liquid nitrogen in a slightly stretched position and stored at −80 °C. Serial cross sections (8-μm thick) were cut in a cryostat microtome set a −20 °C (Thermo Scientific) and incubated with primary antibody against rabbit ClC-1 diluted 1:100 (MyBiosource) in PBS-gelatin for 1 h. After washing with PBS-gelatin, the sections were incubated with 488-donkey anti rabbit IgG (Invitrogen), diluted 1:1000 in PBS-gelatin for 1 h, then washed with PBS-NaCl (300 mM). Sections were examined using Olympus CX41 microscope 53 .

Recordings of Resting Chloride and Potassium Conductances and Excitability Parameters in skeletal muscle of SOD1 G93A and in muscle specific MLC/SOD1 G93A transgenic mice measured in current clamp mode by the two-Intracellular Microelectrodes technique. Extensor Digitorum
Longus (EDL) and Soleus (Sol) muscles were fixed by tendons to a glass rod immersed in normal (NP) or chloride-free physiological solution maintained at 30 °C and perfused with 95% O2/5% CO2 23 . The NP solution contained (in mM): NaCl 148, KCl 4.5, CaCl2 2.0, MgCl2 1.0, NaHCO3 12.0, NaH2PO4 0.44, glucose 5.5, and pH 7.2. The chloride-free solution was prepared by equimolar substitution of methylsulfate salts for NaCl and KCl and nitrate salts for CaCl2 and MgCl2. The cable parameters of myofiber sarcolemma were determined from the electrotonic potentials elicited by square wave hyperpolarizing current pulse of 100-ms duration, using two intracellular microelectrodes in current-clamp mode, as previously described 30,31,54 . The membrane conductance is calculated from the values of input resistance, space constants and time constant and assuming a myoplasmic resistivity of 125 Ω cm. The mean chloride conductance (gCl) is calculated as the mean total membrane conductance (gm) measured in NP solution minus the mean potassium conductance gK measured in chloride-free solution. Sarcolemma excitability parameters were determined by applying 100ms-long depolarizing current pulses of increasing amplitude to elicit first a single action potential (AP) then a train with the maximal number of APs. The membrane potential was held at −80 mV between test pulses. A single and a train of action potentials were generated by increasing current intensity in the same fiber and the excitability characteristics were measured (amplitude of action potential, threshold current, latency of action potential and N spikes) 31,55 . Measurement of KATP channel activity in skeletal muscle of SOD1 G93A and in muscle specific MLC/SOD1 G93A transgenic mice using whole-cell patch-clamp technique. Patch-clamp experiments were performed on freshly enzymatically dissociated FDB muscle fibers in inside-out configurations using the standard patch-clamp technique. The ATP-sensitive potassium (KATP) channel currents were recorded immediately after excision during voltage steps from 0 to −60 mV with 150 mM KCl on both sides of the membrane patches in the absence (control) or the presence of ATP in the muscle bath 24 . The currents were recorded at a 1-kHz sampling rate (filter = 0.2 kHz) using an Axopatch-1D amplifier equipped with a CV-4 head stage (Axon Instruments, Union City, CA). Macropatches having an average pipette area of 11.3 ± 1 μm 2 were used to measure the mean KATP currents, which were calculated by subtracting the baseline level from the open-channel level of each current trace and then digitally averaging all generated files using CLAMPFIT (Axon Instruments). The baseline level for the KATP current was measured in the presence of internal ATP (5 × 10−3 M). Current amplitude was measured using CLAMPFIT.
Fluorescence Measurements of Resting Intracellular Ca 2+ Concentration in skeletal muscle of SOD1 G93A and in muscle specific MLC/SOD1 G93A transgenic mice using FURA-2 imaging technique. Fluorescence measurements were performed on small bundles of five to ten fibers lengthwise dissected from mice EDL muscles, as described elsewhere 56 . The muscle fibers were incubated with the fluorescent calcium probe FURA-2 for 45-60 min at 22 °C in physiological solution containing 5 µM of the acetoxymethyl ester (AM) form of the dye mixed to 10% (v/v) Pluronic F-127 (Molecular Probes, Leiden, The Netherlands). A QuantiCell 900 integrated imaging system (VisiTech International Ltd) was used to acquire pairs of background-subtracted www.nature.com/scientificreports www.nature.com/scientificreports/ images of the FURA-2 fluorescence emission (510 nm) excited at 340 nm and 380 nm. The equation used to transform fluorescence ratio in [Ca 2+ ]i values was [Ca 2+ ]i = (R-Rmin)/(Rmax-R) •K D •β, where R is the ratio of fluorescence excited at 340 nm to that excited at 380 nm; K D = 145 nM; β, Rmin and Rmax were determined in situ in ionomycin-permeabilized muscle fibers. The cytosolic calcium response of muscle fibers to in vitro application of 40 mM caffeine was measured to evaluate calcium release from the sarcoplasmic reticulum. statistical analysis and Multivariate statistical analysis: pCA and LDA algorithms. All data are expressed as the mean ± SEM. Statistical analysis for direct comparison between two groups of data means was performed using unpaired Student's t-test, while multiple statistical comparison between groups was performed using ANOVA test, with Bonferroni's t test post hoc correction, allowing a better evaluation of intraand inter-group variability (using Prism 7.0 Software, GraphPad). For the analysis of correlations between two parameters, Python 3.6, Scipy and Seaborn packages have been used and Pearson correlation coefficient and the 2-tailed p-value for testing non-correlation have been calculated. Multivariate analysis of gene expression dataset has been performed following the PCA-LDA analysis 57,58 . For each animal 27 genes have been analyzed as normalized values ratio with house-keeping gene β-Actin (HK: Actb; target genes: ClC1, PKC Theta, PKC Alpha, Ryr1, Serca1, Serca2, CN, Bk, Kir6.2, Sur1, Surb2a, Sur2b, Nav1.4, Nmdar1, Ampar2, Murf1, Irisin, Notch1, TauT, Ampk, Ngf, Mstn, Achr1, Mhc2b, Mhc2x, Mhc2a, Mhc1). The starting data-frame consisted of a matrix with 30 rows (animals) and 27 columns (gene expression levels) with a total of 810 points. Because PCA and LDA algorithms suffer of the presence of big outliers in the dataset, first step of the pipeline has been their handling. So, the two biggest outliers onto genes that represents the major part of the variance in the system of data (ESD method, p < 0.001) have been removed and the resulting empty cells of table have been filled with the mean values of the mRNA's level of the groups to which they belong. The second step has been to standardize the gene expression dataset scaling the mean to 0 and the variance to unitary scale, this is done for each column. After that a PCA object has been trained by the use of the Scikit-learn class sklearn.decomposition.PCA setting the number of principal component computed up to 27 (the minor dimension of the dataset) and the svd-solver parameter to 'full' , because the dataset has less than 500 × 500 dimension. In this way, the object PCA has been trained, taking as input, the consolidated and standardized dataset, computing than a linear dimensionality reduction using the Singular Value Decomposition solver. The subsequent step in the PCA-LDA pipeline was to perform 3 Fisher Linear Discriminant Analysis (LDA) of PC scores. The class sklearn.discriminant_analysis has been used to build a prototype classifier with a linear decision boundary using the Bayes' rule, starting from a matrix n x m and a vector containing group identifiers. The 3 LDA further reduced dataset's dimension starting the first three PC scores, to an single Linear Discriminant (LD) dimension, in a supervised manner, maximizing the variance between and minimizing the variance within classes. From each LDA a 1 × 3 vector with the LD coefficients of the 3 PC have bene calculated. In the first LDA we obliged the classifier to find the LD that better separate MLC/ SOD1 G93A mice from WT. In the second LDA the LD that better separate SOD1 G93A mice from the WT and in the third LDA the LD that better segregate MLC/SOD1 G93A mice from the SOD1 G93A ones. Finally, 3 vectors of LDA coefficients and a matrix containing the PCA loading weights of each genes have been obtained. Thus, the three discriminant lists of genes that separate the classes have been obtained multiplying the PCA loadings matrix with each transposed vector of LDA coefficients. Each multiplication gave as result, a vector of 27 PCA-LDA discriminant weights of genes analyzed. The absolute values of each element in these lists have been ordered in growing manner so the most important discriminant genes appeared at the bottom of the three lists of discriminant genes.