Transcriptome Profiling of Spinal Muscular Atrophy Motor Neurons Derived from Mouse Embryonic Stem Cells

Proximal spinal muscular atrophy (SMA) is an early onset, autosomal recessive motor neuron disease caused by loss of or mutation in SMN1 (survival motor neuron 1). Despite understanding the genetic basis underlying this disease, it is still not known why motor neurons (MNs) are selectively affected by the loss of the ubiquitously expressed SMN protein. Using a mouse embryonic stem cell (mESC) model for severe SMA, the RNA transcript profiles (transcriptomes) between control and severe SMA (SMN2+/+;mSmn−/−) mESC-derived MNs were compared in this study using massively parallel RNA sequencing (RNA-Seq). The MN differentiation efficiencies between control and severe SMA mESCs were similar. RNA-Seq analysis identified 3,094 upregulated and 6,964 downregulated transcripts in SMA mESC-derived MNs when compared against control cells. Pathway and network analysis of the differentially expressed RNA transcripts showed that pluripotency and cell proliferation transcripts were significantly increased in SMA MNs while transcripts related to neuronal development and activity were reduced. The differential expression of selected transcripts such as Crabp1, Crabp2 and Nkx2.2 was validated in a second mESC model for SMA as well as in the spinal cords of low copy SMN2 severe SMA mice. Furthermore, the levels of these selected transcripts were restored in high copy SMN2 rescue mouse spinal cords when compared against low copy SMN2 severe SMA mice. These findings suggest that SMN deficiency affects processes critical for normal development and maintenance of MNs.


Introduction
Spinal muscular atrophy (SMA) is an autosomal recessive, earlyonset neurodegenerative disorder characterized by the degeneration of a-motor neurons (MNs) in the anterior horn of the spinal cord which leads to progressive muscle weakness and atrophy [1]. SMA is a leading genetic cause of infant death worldwide with 1 in 5000-10,000 children born with the disease [2,3] and a carrier frequency of 1:25-50 [4][5][6][7]. SMA results from the loss or mutation of the SMN1 (survival motor neuron 1) gene on chromosome 5q13 [8]. There is an inverted duplication of SMN1 in humans called SMN2 [9,10]. The duplication of SMN1 only occurs in humans. Within SMN2, there is a C-to-T transition in an exonic splice enhancer of exon 7 that results in the vast majority (about 80-90%) of SMN2 mRNAs to lack exon 7 (SMND7). SMND7 is not fully functional and prone to degradation [11,12]. SMN2 can, however, provide some fully functional SMN protein. The copy number of SMN2 modifies the severity of SMA phenotype in humans [13][14][15][16][17][18][19][20]. Although the genetics underlying SMA are known, the mechanisms leading to the disease are poorly understood.
SMN is a ubiquitously expressed protein that facilitates the assembly of ribonucleoproteins (RNPs) [21]. The assembly of Utype small nuclear RNPs (snRNPs) is developmentally regulated in many cell types [22]. Knockdown of snRNP assembly in zebrafish results in degeneration of motor neuron axons [23]. snRNP assembly is defective in SMN-deficient SMA cells [24]. However, this function of SMN in U-type snRNP assembly is required by all cells. So why are MNs primarily affected in SMA given this ubiquitous function? Gabanella et al. show that snRNP assembly is defective in tissues from mouse models for SMA and that the extent of reduced snRNP assembly correlates with phenotypic severity of these SMA mice [22]. Furthermore, snRNP assembly is more markedly affected by Smn deficiency in SMA mouse neural tissues than in other tissues like the kidney. This enhanced sensitivity of neurons to deficits in snRNP assembly may explain why MNs are primarily affected in SMA.
In addition to its role in snRNP biogenesis, SMN may also have functions which are unique to MNs. SMN associates with cytoskeletal elements in MN dendrites and axons within the rat spinal cord [25]. Knockdown of Smn expression in zebrafish embryos using morpholino antisense oligonucleotides leads to reduced axonal outgrowth and axonal pathfinding defects [26]. Ectopic overexpression of wild-type SMN in these Smn-knocked down zebrafish embryos restores normal axonal growth and pathfinding. Moreover, axonal defects in Smn-knocked down zebrafish embryos are corrected by overexpression of mutant SMNs which are incapable of snRNP assembly [27]. The core protein components of snRNPs, the Sm proteins, do not colocalize with SMN in neuronal processes [28]. Taken together, these observations suggest that SMN may have a unique function in neurons that is independent of snRNP biogenesis. Several studies have suggested that SMN may play a role in axonal trafficking. SMN has been shown to interact with the heterogeneous nuclear ribonucleoprotein-R (hnRNP-R) [29]. hnRNP-R binds to the 39untranslated region (39-UTR) of b-actin mRNA [30]. Transport of b-actin mRNA and protein into growth cones has been shown to be essential in developing neurites [31]. Reduced expression of Smn in mouse MNs results in reduced axonal outgrowth and trafficking of b-actin mRNA along axons [30]. Conditional knockout of b-actin in mouse MNs, however, does not affect the morphologies of MN neurites suggesting that b-actin mRNA transport along axons may not be essential for MN axon development [32]. Although neuron-specific roles of SMN have been proposed, further studies are needed to determine how these neuron-specific roles are involved in the pathogenesis of SMA.
Mouse models have been instrumental in studying the function of SMN and how this function is disrupted in SMA. Mice, like all animals except humans, have only one SMN gene (mSmn) that is equivalent to SMN1. Homozygous disruption of mSmn in mice is embryonically lethal [33]. Transgenic overexpression of SMN2 rescues the embryonic lethality of mSmn deficiency [34,35]. mSmn nullizygous mice that also harbor few copies of the SMN2 transgene (i.e. 1 or 2) develop a severe motor phenotype resembling SMA and die within 7 days after birth [34,35]. Increasing the SMN2 copy number in these mSmn nullizygous mice improves the survival and phenotype of these SMA mice; in fact, expression of 8-16 copies of SMN2 fully rescues the SMA phenotype in these mice [34,36]. Patients who have been identified genetically as SMA-i.e. loss of SMN1-are phenotypically normal when they carry at least 5 copies of SMN2 [16]. Thus, SMN2 expression modifies the phenotypic severity of SMA in mice as well as in man and makes SMN2 a target for the development of SMA therapeutics.
The low copy SMN2 SMA mouse phenotypically resembles type I SMA in humans [34]. The short lifespan as well as the low frequency of pups that survive past birth limit their use for mechanistic studies; therefore, an in vitro model would be useful for such studies. Murine embryonic stem cells (mESCs) are able to differentiate into spinal neural progenitor cells and then into MNs through exposure to the morphogens retinoic acid (RA) and Sonic hedgehog (Shh) [37]. Motor neurons differentiated from mESCs were found to generate action potentials and developed axons and synapses when co-cultured with muscle cells [38]. mESC lines have been established for low copy SMN2 severe SMA mice also harboring a MN-specific reporter construct (HB9:eGFP) [39]. When these SMA mESCs are directed to differentiate into MNs, they start dying after the differentiation process [39]. MNs derived from SMA mESCs can, therefore, potentially provide important insights into the pathogenesis of SMA.
In this study, we will use cultured MNs derived from SMA mESCs to determine how reduced levels of the ubiquitously expressed protein SMN result in selective MN death in SMA. Previous studies have used cDNA microarrays to identify differentially expressed mRNAs in SMA mouse whole spinal cords and in primary MN cultures [40][41][42]. Microarrays can only identify known RNA transcripts which limits their utility for comprehensively characterizing transcriptomes. Massively parallel RNA sequencing, commonly known as RNA-Seq [43], is a recently developed deep-sequencing technology used for transcriptome profiling [44]. RNA-Seq directly reads the sequences of the cDNA pool which results in a very low background signal as compared to the indirect method of measuring hybridization intensity used in microarray analysis. Since RNA-Seq directly reads cDNA sequences, novel transcripts and isoforms can be identified. In this study, we use RNA-Seq to annotate and compare the transcriptome profile of MNs derived from severe SMA mESCs with those derived from normal mESCs. Analysis of over-represented biological pathways and networks revealed that SMA mESC-derived MNs have increased expression of RNA transcripts related to pluripotency and reduced expression of neuronal development and function RNA transcripts. This study provides new insights into the molecular consequences of SMN deficiency in MNs and identifies novel targets for the development of neuroprotective therapeutics.

Ethics Statement
All animal experiments were conducted in accordance with the protocols described in the National Institutes of Health Guide for the Care and Use of Animals and were approved by the Nemours Biomedical Research Institutional Animal Care and Use Committee.
Spinal cords were collected at 2 time points: embryonic day 13.5 (e13.5) and postnatal day 3 (PND03; with the day of birth being defined as PND01). For collecting e13.5 samples, timedpregnant dams (as assessed by the presence of a vaginal plug) were euthanized at e1360.5 and the spinal cords were rapidly dissected from the embryos, snap-frozen and stored at 280uC until RNA isolation. Additional tissues were harvested from each embryo for genotyping. For postnatal samples, pups were euthanized and the spinal cords were rapidly dissected from the pups, snap-frozen and stored at 280uC until RNA isolation. Tail biopsies were also taken for genotyping. All embryonic and postnatal samples were genotyped as described previously [34,48].

Fluorescence-Activated Cell Sorting
Dissociated ES cells were filtered through a cell strainer, suspended in DMEM with 20 mg/mL of propidium iodide (PI) and sorted using a 488 nm laser attached to either a MoFlo cell sorter (Coulter; Kimmel Cancer Center, Thomas Jefferson University) or the FACSAria II sorter (BD Biosciences; Center for Translational Cancer Research, University of Delaware). Cells were passed through a 100 mm nozzle tip at a speed of approximately 12,500 events per second. Embryonic stem cells from nonGFP-expressing, wild-type mice described previously [46] were used as a negative control to set the cutoff for background fluorescence. Cells that were GFP positive, PI negative were collected in a 5 mL round bottom test tube (BD Biosciences, Franklin Lakes, NJ, USA) with PBS supplemented with 60% FBS. Approximately a million cells were collected per sample. Sorted cells were centrifuged, PBS was removed and pellets were snap frozen for RNA isolation.

RNA Isolation
Total RNA was isolated from cells using the RNeasy mini kit (QIAGEN, Germantown, MD, USA) with an additional DNase I (QIAGEN) digestion step to remove any genomic DNA contamination. The concentration of the purified RNA was determined by a ND-1000 NanoDrop the spectrophotometer (NanoDrop Technologies, Wilmington, DE). RNA integrity was assessed by the Agilent Technologies 2100 Bioanalyzer.

RNA-Seq
1 mg of RNA from each sample (n = 3/genotype) was collected and sent to the Sequencing and Genotyping Center at the Delaware Biotechnology Institute in the University of Delaware for RNA-Seq library preparation and sequencing. The cDNA library was prepared using the TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA) according to the manufacturer's directions. The samples were then clustered and sequenced on an Illumina HiSeq 2500. Deep sequencing was performed on triplicates for each cell line (total six samples) for a 50 cycle single end run. The raw sequence data have been deposited in the NCBI Gene Expression Omnibus (GEO) [49] under accession number GSE56284.

RNA-Seq Data Analysis
RNA-Seq reads were assessed for quality control with FastQC (version 0.10.1; Babraham Bioinformatics, Cambridge, UK). Adapter sequences and poly-A tails were removed from sequencing reads with CutAdapt (MIT, Cambridge, MA) [50] and were confirmed by FastQC. Reads were mapped to a reference mouse transcriptome and genome (NCBI m37) using TopHat (version 2.04; [51]). The reference genome and annotation was obtained from Ensembl (http://www.ensembl.org). Transcripts were assembled and transcript abundances were measured as fragments per kilobase of exon per million fragments mapped (FPKM) by Cufflinks [52]. Cuffdiff [53] was then used to determine differential expression. Differential expressed genes were plotted using the R package CummeRbund [53].

Gene Pathways and Networks Analysis
Identification of biological pathways and networks affected by Smn deficiency in mESC-derived MNs was completed using Ingenuity Pathway Analysis (IPA build version 261899; Ingenuity Systems, Inc., Redwood City, CA). The list of differentially expressed transcripts generated from Cuffdiff was divided into upregulated and downregulated lists of transcripts and were imported into IPA. Biological function and canonical pathways were determined to be over-represented using the Fisher exact test with a false discovery rate (FDR) correction (p#0.05). A set of interacting biological networks was generated in IPA from the upregulated and downregulated transcript lists. A network score was computed based on the likelihood that the focus transcripts specifically belong to this network as opposed to being in the network by chance. Based on other published work [54], only those networks with scores greater than or equal to 15 were considered relevant for further analysis.
The Upstream Regulator Analysis (URA) algorithm of IPA allows for the identification of upstream regulators-transcription factors, signaling molecules and drugs-which can affect gene expression [55]. URA was completed on upregulated and downregulated transcripts and list of putative upstream regulators was provided. Each upstream regulator was assigned an activation z-score; upstream regulators were considered as being activated or inhibited if their activation z-scores were either greater than or equal to 2.0 for activated regulators or less than or equal to 22.0 for inhibited regulators.
cDNA Synthesis and Quantitative Reverse Transcriptase-Polymerase Chain Reaction (qRT-PCR) cDNA was prepared from 0.5-1 mg total RNA using the iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA) as per manufacturer's instructions. The cDNA-diluted 200-400 foldwas amplified via RT-PCR using the QuantiTect SYBR Green PCR kit (QIAGEN). The target primers listed in Table 1 were synthesized by Integrated DNA Technologies (Coralville, IA) and primers for the reference transcripts b-glucuronidase (Gusb), phosphoglycerate kinase 1 (Pgk1) and ribosomal protein L13a (Rpl13a) were purchased from Real Time Primers LLC (Elkins Park, PA). Quantitative PCR was performed in a 384 well plate on a 7900 HT Fast Real-Time PCR system (Applied Biosystems, Foster City, CA). Each sample was assayed in triplicate. Relative transcript levels were calculated using the comparative C t (2 2DDCt ) method [56] where DC t is defined as the difference between the C t for the target transcript and the C t for the geometric mean of Gusb, Pgk1 and Rpl13a [57] and DDC t is defined as the difference between the DC t for the SMA sample and the DC t for the control sample.

Statistical Analysis
All quantitative data is expressed as mean 6 standard error. For differential expression of RNA transcripts analyzed by RNA-Seq, statistical analysis was completed using Cuffdiff and CummeRbund [53]. Statistical analysis for pathways and networks analysis was completed by IPA. Statistical analyses of quantitative data were completed using SPSS v.22.

Control and SMA mESCs Differentiate into Motor Neurons in vitro
Motor neurons (MNs) were generated from HB9:eGFPexpressing control (SMN2 +/+ ;mSmn +/+ , Hb9) and SMA (SMN2 +/ + ;mSmn 2/2 , A2) mESCs using a combination of growth factors and morphogens [47] as outlined in Figure 1. We first examined the effect of MN differentiation of the levels total SMN protein levels (mSmn and the SMN2 transgene) in Hb9 and A2 cells. Total SMN protein levels were indeed reduced in undifferentiated as well as in differentiated A2 cells relative to Hb9 cells ( Figure 2). The total SMN protein levels were similar between undifferentiated and differentiated Hb9 as well as for A2 cells; the differentiation conditions thus do not affect total SMN protein expression.
We next determined the effect of SMN deficiency on the differentiation potential of mESCs into MNs. FACS analysis showed the differentiation efficiency of viable control Hb9 cells ( Figure 3B; 17.766.1%) was similar to that observed in viable A2 SMA cells ( Figure 3C; 16.363.2%). Reduced expression of SMN similar to that observed in SMA does not affect the ability of mESCs to differentiate into MNs.
To show that the differentiated cells from Hb9 and A2 mESCs were morphologically MNs, these mESCs were directed to differentiate into MNs and then plated in the presence of neurotrophic factors to induce axonal growth and arborization ( Figure 1). The GFP+ Hb9 (control in Figure 4) and A2 (SMA in Figure 4) cells were characterized using various neuronal markers to test their MN identity. The presence of b-III tubulin (Tuj1), a neuron-specific tubulin isoform, as well as post-mitotic neuronal marker NeuN (also known as Fox-3) confirmed that the GFP+ cells were indeed neurons ( Figure 4A). As expected, differentiated Hb9 and A2 cells expressed MN specific markers Hb-9 and Islet-1 ( Figure 4B) [58,59]. One noticeable difference was the number of GFP+ neuron-like cells present on coverslips after plating and growing for 3 days. Coverslips with differentiated Hb9 cells had an abundance of GFP+ cells whereas there were very few neuron-like GFP+ cells on coverslips from A2 SMA cells. This reduced viability of SMA mESC-derived MNs was consistent with a previous report [39]. Thus, both the control and SMA mESCs can be directed to differentiate into MNs and SMNdeficient SMA MNs are less viable than control MNs derived from mESCs.   Analysis of RNA-Seq Data from Control and SMA mESC-Derived MNs Whole transcriptome differential expression between FACScollected Hb9 and A2 MNs was completed using massively parallel RNA sequencing (RNA-Seq). The total number of reads produced from each sample (n = 3/genotype) was around 100 million reads per sample. The quality of the trimmed reads-as assessed with FastQC-was very high in that the sequence quality had a Phred score of around 39 and an average base Phred quality score was around 32. As a point of reference, a Phred quality score of 30 corresponds to a base call accuracy of 99.9%.
The trimmed sequencing reads were then mapped to the mouse reference genome ( Figure 5A). The total number of mapped reads per sample ranged between 29,763,880 and 44,570,352 ( Table 2); between 69% and 76% of reads were uniquely mapped to the reference genome. The difference in the average number of reads between Hb9 and A2 MNs samples was not statistically significant (p = 0.22). Transcripts were assembled and transcript abundance measured in FPKM; the expression levels for each transcript were then compared between Hb9 and A2 MNs. The distributions of transcript abundances-as measured in fragments per kilobase of exon per million fragments mapped (FPKM)were different between A2 SMA and Hb9 control MNs ( Figure 5B). Overall, there were 41,945 expressed transcripts in both Hb9 and A2 MNs of which 10,058 transcripts were found to be differentially expressed in Hb9 and A2 MNs (p,0.05). Out of the 10,058 statistically significant differentially expressed transcripts, 3094 were upregulated in A2 SMA samples compared to Hb9 control and 6964 were downregulated ( Figure 5C). A list of differentially expressed RNA transcripts is provided in Table S1.
Since the A2 SMA cells are homozygous for the targeted deletion of mSmn (Smn1), the levels of Smn1 mRNA is expected to be lower in A2 SMA MNs. Indeed, the level of Smn1 mRNA was reduced by 4.7-fold in A2 SMA MNs when compared to Hb9 controls MNs (Table S1). The mRNA transcript data (RNA-Seq) and the protein results (Smn immunoblot in Figure 2) were consistent. Interestingly, Hb9 (Mnx1) transcripts levels were not significantly different between A2 and Hb9 MNs. This observation was supported by the FACS analysis in Figure 3.

Pathway and Network Analysis of Differentially Expressed Transcripts in SMA ESC-derived MNs
Ingenuity Pathway Analysis (IPA) was used to determine the biological relevance of the differentially expressed transcripts between A2 SMA and Hb9 control mESC-derived MNs. IPA uses the Ingenuity Knowledge Base, which is a manually curated repository of literature-based biological and pharmacological data (www.ingenuity.com) [55]. 77 biological functions were overrepresented (with a Fisher's exact test p-value cutoff of 0.05) using the list of transcripts upregulated in A2 SMA MNs, the 7 highest scoring entries are shown in Figure 6A. Many of these highest scoring entries are related to developmental functions. In the downregulated transcripts, there were 62 overrepresented biological functions ( Figure 6B) with many of these highest scoring functions involving nervous system development and disease.
We also used IPA to identify canonical pathways that were overrepresented in the lists of transcripts either upregulated or downregulated in A2 SMA MNs. Most of the highest scoring overrepresented canonical pathways (6 out of 14) from the upregulated transcripts were related to ESC pluripotency ( Figure 6C with the top canonical pathway, Transcriptional Regulatory Networks in ESCs, shown in detail in Figure 6E). Many of the highest scoring canonical pathways identified from the list of downregulated transcripts (5 out of 19) were involved with nervous system development ( Figure 6D). In fact, the top downregulated canonical pathway is Notch signaling (shown in detail in Figure 6F).
In addition to identify biological functions and canonical pathways altered as a result of SMN deficiency in MNs, IPA can also integrate these functions and pathways to identify functionally interacting gene networks, or interactomes. We analyzed our lists of transcripts upregulated or downregulated in A2 SMA MNs for overrepresented interactomes. IPA network scores greater than or equal to 15 were considered significant [54]. Using these criteria, there were 2 significant interactomes identified from analysis of the upregulated transcripts and 13 from analysis of the downregulated transcripts. The top upregulated interactome (score of 59; Figure 7A) contained many gene products involved in ESC pluripotency and cellular development, growth and proliferation. Gene products with roles in nervous system development and function were present in the top scoring interactome of the downregulated transcripts (score of 60; Figure 7B). Taken together, a systems biology-based comparison of A2 SMA and Hb9 control mESC-derived MNs shows that functions, pathways and gene networks involved in pluripotency are upregulated in SMA MNs while those involved in nervous system development and function are downregulated in SMA MNs.

In Silico Identification of Upstream Regulators in SMA mESC-derived MNs
Using the Upstream Regulator Analysis (URA) component of IPA, we can identify possible upstream regulatory molecules that may be affecting changes in gene expression [55]. A list of such upstream regulators was compiled by URA using the transcripts which were upregulated or downregulated in A2 SMA mESCderived MNs. We divided these upstream regulators into 3 groups: transcriptional regulators, signaling components and drugs. As shown in Figure 8A, 4 transcriptional regulators were activated in SMA MNs while 13 transcriptional regulators were inhibited in SMA MNs. Of these 17 transcriptional regulators identified, the transcript levels of only three-Pou5f1 (Oct-4), Nanog and Ascl1-were significant altered in A2 SMA MNs (3.0-fold increase, 2.5-fold increase and 5.6-fold decrease, respectively).
In addition to transcriptional upstream regulators, the intracellular effects of 18 signaling molecules were identified as being significantly modified in A2 SMA MNs ( Figure 8B). These signaling molecules include growth factors and cytokines like transforming growth factor b (TGFb) and neurotrophin-4 (NTF-4), intracellular second messengers like cyclic AMP (cAMP) as well as nuclear receptor ligands such as dihydrotesterone (DHT) and tretinoin (RA). The effects of all of these signaling molecules except for inhibin bA (INHBA) were attenuated in A2 SMA MNs. The levels of Sonic hedgehog (Shh) mRNA were significantly reduced (2.4-fold) in A2 SMA MNs relative to Hb9 control MNs.

Validation of Differentially Expressed Transcripts Identified by RNA-Seq
As a first step toward the validation of the RNA-Seq data, we needed to reduce the number of differentially expressed transcripts to a more manageable list. The initial list of differentially expressed transcripts (using a p#0.05 as a criterion) was re-analyzed using more stringent criteria. By filtering these data for differentially expressed transcripts with a p-value less than or equal to 0.01 and a greater than 4-fold change, the list of candidate transcripts was reduced to 286 upregulated transcripts and 814 downregulated genes. The filtering of these data was further refined so as to include only those transcripts with a FPKM .30 ( Figure 5B). With this added refinement, there were 27 transcripts upregulated in A2 SMA MNs and 220 downregulated transcripts.
To validate the results from the analysis of the RNA-Seq data, we measured changes in the levels of selected differentially expressed transcripts between Hb9 and A2 MNs by qRT-PCR. We selected Smn1 since this gene is knocked out in SMA A2 cells. In addition, six other transcripts were selected based on their strong changes in transcript levels as shown by RNA-Seq: cellular retinoic acid binding protein 1 (Crabp1), Crabp2, Islet-1 (Isl1), NK2 homeobox 2 (Nkx2.2), phospholipase A2, group 1B (Pla2g1b) and vimentin (Vim). The sample RNAs used for qRT-PCR (n = 3/ genotype) were not the same as those used for RNA-Seq so they represent biological replicates as opposed to technical replicates [60]. The differences in transcript levels between Hb9 and A2 MNs determined by qRT-PCR followed the same trends as those determined by RNA-Seq although the magnitudes of change were generally higher in the RNA-Seq data ( Figure 9A). RNA-Seq is more sensitive than qRT-PCR at detecting changes in transcript levels.
We next determined if the changes in RNA levels observed in these SMA mESC-derived MNs were unique to these specific cells. Control and severe SMA mESCs that do not contain the HB9-eGFP reporter transgene-C4 and E2 cells, respectively- [46] were directed to differentiate into MNs. The extracted total RNAs from C4 and E2 MNs were analyzed by qRT-PCR. As shown in Figure 9B, with the exception of Nkx2.2, the changes in transcript levels between C4 and E2 MNs were similar to those observed between eGFP-labeled Hb9 and A2 MNs. One possible explanation for this discrepancy in changes in Nkx2.2 levels is that the MNs were not separated from other cell types after differentiation of C4 and E2 mESCs.
Using two-dimensional differential in-gel electrophoresis (2D-DIGE) and mass spectrometry, we previously identified a set of proteins that are differentially expressed between control and severe SMA mESC-derived MNs [46]. We compared this list of differentially expressed proteins with our list of differentially expressed mRNA transcripts. Of the 12 proteins identified, lactate dehydrogenase B (Ldhb) and aldehyde dehydrogenase (Aldh5a1) showed corresponding changes in mRNA levels between control and SMA mESC-derived MNs ( Table 3). 5 proteins that showed significant difference in protein levels between control and SMA mESC-derived MNs-brain creatine kinase (Ckb), tropomyosin 3 (Tpm3), ubiquitin C-terminal hydroxylase L1 (Uchl1), 14-3-3c (Ywhag) and heat shock protein 90b (Hsp90b1)-did not exhibit any corresponding changes in mRNA levels. While some genes have similar changes in mRNA and protein expression in SMA MNs relative to control MNs, there are other genes which are differentially expressed at the protein level but not at the mRNA level. These observations suggest that differential gene regulation occurs at many levels-i.e. transcriptional vs. post-transcriptional-between SMA and control MNs.

Differential Expression of Validated Transcripts in SMA Mice
The levels of Smn1, Crabp1, Crabp2, Isl1, Nkx2.2, Pla2g1b and Vim transcripts were examined in total RNA samples from control (SMN2 +/+ ;mSmn +/+ ) and severe SMA (SMN2 +/+ ;mSmn 2/2 ) mouse spinal cords in order to determine if the changes observed in mESCderived MNs could also be observed in vivo. Mouse embryos of similar genotypes were used to generate the mESCs used in this study. Spinal cord total RNAs (n = 3/genotype) were extracted from PND03 mice; severe SMA mice at this time point begin to show signs of motor dysfunction [34]. Similar to SMA mESC-derived MNs in culture, Smn1, Crabp1, Crabp2 and Nkx2.2 transcript levels were reduced while Pla2g1b levels were increased in SMA spinal cords ( Figure 10A). Surprisingly, Isl1 and Vim mRNA levels were elevated in SMA spinal cords at PND03 even though these transcripts were reduced in SMA MNs. The samples isolated from SMA mouse spinal cords contain RNAs from many different types of neurons aside from MNs as well as other cell types such as astrocytes and oligodendrocytes. This sample heterogeneity could explain the discrepancies observed between mESC-derived SMA MNs and SMA spinal cords.
Increasing SMN2 copy numbers can improve the phenotype and survival of severe SMA mice [34,36]. In fact, SMN2 transgenic SMA mice with 8 (SMN2(566) +/2 ;mSmn 2/2 ) or 16 copies (SMN2(566) +/+ ;mSmn 2/2 ) of the transgene display no motor phenotype [34]; in other words, the SMA phenotype is rescued. When comparing relative changes in Smn1, Crabp1, Crabp2, Isl1, Nkx2.2, Pla2g1b and Vim transcript levels in lowcopy SMN2 SMA mice with those observed in high-copy SMN2 rescue mice at PND03, the relative changes in each of these transcripts-except for Smn1-were increased in high-copy rescue mice as opposed to relative changes in low-copy SMN2 SMA mice ( Figure 10A). Smn1 mRNA levels were still reduced in high-copy SMN2 rescue mice since these mice are still nullizygous for mSmn.
We also wanted to determine if the changes observed in severe SMA mouse spinal cords are regulated developmentally. The levels of Smn1, Crabp1, Crabp2, Isl1, Nkx2.2, Pla2g1b and Vim transcripts were measured in severe SMA mice at embryonic day 13.5 (e13.5) relative to age-matched control littermates. There are no overt changes in physical characteristics or in motor neuron growth in SMA mice at this developmental age [45]. Smn1 mRNA levels were marked reduced in the embryonic SMA mouse spinal cord as they were at PND03 ( Figure 10B). Some transcripts, like Isl1, Pla2g1b and Vim, showed reduced mRNA levels at e13.5 even though the levels of these transcripts were increased at PND03. The reductions in Crabp2 and Nkx2.2 transcripts were more pronounced at e13.5 than at PND03. These results suggest that the changes in RNA expression of some transcripts may be developmentally regulated in SMA.

Discussion
In this study, we compared the RNA expression profiles (transcriptomes) between SMA and normal MNs in order to identify the molecular consequences of SMN deficiency in these cells. This study is the first to examine the transcriptome profiles of isolated MNs derived from severe SMA and normal mESCs. Differential gene expression studies using cDNA microarrays have also been previously completed using primary MN cultures from mSmn +/2 mice-a supposed model for mild SMA [61]-and whole spinal cords from the SMA mouse models [40][41][42]. As mentioned earlier, analysis of transcriptome profiling using RNA-Seq is advantageous because of the high signal-to-noise ratio, independence from hybridization efficiency between probes and targets and identification of novel RNA transcripts. In addition to mouse models for SMA, differential transcriptome expression analyses have been completed on other animal models of Smn  deficiency [62,63]. Recently, Zhang et al. [64] used RNA-Seq to identify differentially expressed mRNAs in PND01 SMND7 SMA mouse MNs isolated by laser capture microdissection (LCM). Our study examines the transcriptome in purified, cultured MNs derived from control and severe SMA mESCs. Using this model, we can examine changes in RNA expression in a nearly homogeneous population of cells which is not possible with LCM.An advantage of using MNs derived from mESCs to compare the transcript profiles of SMA MNs to normal MNs is that large numbers of cells can be directed to differentiate into MNs. This scalability can overcome two obstacles of using primary MN cultures from SMA mouse spinal cords, low yield of MNs from SMA mouse embryos and the fragility of these affected MNs.
Numerous studies in mice and zebrafish have shown that SMA MNs are affected in a cell autonomous manner [26,[65][66][67] which would suggest that examining the transcriptome profiles of isolated SMA MNs could provide relevant information on the pathobiology of this disease. In isolated mESC-derived SMA MNs, the transcriptomes will not show the effects of target (i.e. skeletal muscle) innervation and trophic signaling from skeletal muscle, sensory neurons or glia. Future experiments will compare the transcriptomes of mESC-derived SMA MNs grown in isolation with those grown in the presence of myotubes and/or glial cells so as to determine those gene regulatory events which are intrinsic to SMA MNs (i.e. those grown in isolation) and those which are dependent on environmental cues (i.e. those grown in the presence of target or support cells).
MNs are the primary cells affected by reduced SMN expression in SMA. Ectopic overexpression of SMN in the neurons of severe SMA mice rescues the primary disease phenotype in these mice while transgenic overexpression of SMN in mature skeletal muscle does not improve the SMA phenotype [68]. Conditional expression of SMN in the developing MNs of SMA mice-using either the Hb9 or Olig2 promoters as drivers-significantly ameliorates the SMA phenotype [65][66][67]. Martinez et al. [66] also show that conditional expression of SMN in SMA skeletal muscle may help grow and maintain muscle independent of MNs. Increasing SMN expression outside of the nervous system with either splice-switching oligonucleotides [69] or adeno-associated virus (AAV) vectors [70,71] markedly improves the phenotype and survival of SMA mice. These studies suggest that comparative analysis of SMA MN transcriptomes from these models may provide limited insight into the pathobiology of SMA; however, it is appropriate to examine the transcript profiles of isolated SMA MNs since they are affected in a cell autonomous fashion [26,[65][66][67].
Many of the biological pathways and networks that were overrepresented in those transcripts upregulated in A2 SMA MNs involved ESC pluripotency ( Figures 6A and 6C). The transcription factors Nanog, Pou5f1 (Oct4), and Sox2 are considered to be hallmarks of ESC pluripotency [78]. mRNA transcripts for all three of these factors were upregulated in SMA mESC-derived MNs ( Figure 6E). UPA of the differentially expressed transcripts revealed that these three pluripotency transcription factors were activated in A2 SMA mESC-derived MNs ( Figure 8A). Several gene products work with these three transcription factors to regulate pluripotency in ESCs. Klf2 regulates the expression of Sox2 [79]. Klf2 transcript levels were increased in SMA mESCderived MNs by 2.3-fold (Table S1). Zic3-whose transcript levels were increased 3.1-fold in SMA mESC-derived MNs-is directly regulated by all three transcription factors [80]. Zscan10 (also known as Zfp206), whose mRNA levels are elevated by 2.5fold in SMA mESC-derived MNs (Table S1), helps maintain pluripotency by jointly functioning with Sox2 and Oct4 [81]. In SMA mESC-derived MNs, the pluripotency marker Dppa5 (also  Table S1). The levels of Dppa5 mRNA are significantly upregulated-as shown with microarray analysis-in the presymptomatic severe SMA mouse spinal cord [42]. Two genes that are also implicated in the maintenance of pluripotency-Mybl2 and Zfp42 (Rex1) [83,84]-were upregulated by 2.2-and 2.8-fold, respectively, in SMA mESC-derived MNs (Table S1).
Neuronal development and activity functions are impaired in A2 SMA mESC-derived MNs given that these biological pathways and networks were overrepresented in the downregulated RNA transcripts (Figures 6B and 6D). Vimentin (Vim) is an intermediate filament protein that is transiently expressed by all neuronal precursors. It may have a role in neurite outgrowth in developing hippocampal neurons [85]. Nkx2.2 is homeobox transcription factor which is activated downstream of Shh; it has a primary role in ventral neuronal patterning and in generating interneurons [86]. Islet1 (Isl1) and Olig2 are transcription factors expressed during and are necessary for MN development [37]. Vim, Nkx2.2 and Isl1 mRNA levels were markedly reduced in mESC-derived SMA MNs ( Figure 9) and in e13.5 SMA mouse spinal cords ( Figure 10); Olig2 mRNA levels were reduced by 5.5-fold in A2 SMA mESC-derived MNs (Table S1). Using microarray analysis, the expression of neurexin2a (nrxn2a) is downregulated in zebrafish embryos with reduced smn expression [63]. Nrxn2a transcript levels are also reduced in the spinal cord of severe SMA mice [63]. Nrxn2 deficiency has been implicated in altered neuronal excitability. We found that Nrxn2a mRNA levels were reduced by 2.9-fold in SMA mESC-derived MNs ( Table  S1).
The Notch pathway, via cell-cell interactions, stimulates neural precursors cells to differentiate into neurons and glial cells [87]. Notch pathway signaling was the top canonical pathway overrepresented in the list of transcripts downregulated in A2 SMA mESC-derived MNs ( Figure 6D). Many of the components of this pathway including Notch as well as its extracellular ligands Jagged and Delta were reduced in A2 SMA mESC-derived MNs ( Figure 6F). Neurogenin-3 (Neurog3) mRNA levels were reduced by 1.3-fold in SMA mESC-derived MNs (Table S1) and Neurog3 was identified as one of the inhibited transcription factor activities in the UPA ( Figure 8A). These results are surprising given that Notch protein immunolocalization is increased on the ventral horn motor neurons of SMND7 SMA mice at PND11 [88]. Jagged1 and Delta1 protein immunolocalization is also increased on adjoining astrocytes in these mice. The differences between our observations and those of Carabello-Miralles et al. [88] could be the result of different SMA models being used (low copy SMN2 severe SMA vs. SMND7 SMA) or types of gene regulation being examined (transcription vs. translation). Cell-cell interactions between neurons and glia that regulate the expression of Notch signaling ligands and receptors may not be apparent in our model system (purified MNs in culture); however, both studies document reduced Neurog3 expression and activity in SMA MNs ( [88] and this work).
RA is known for driving cellular differentiation and has a neuralizing effect during neuronal development. UPA of the differentially expressed transcripts in A2 SMA mESC-derived MNs showed that all trans-RA (tretinoin) signaling was inhibited in SMA MNs ( Figure 8B). Transcripts whose protein products are involved in RA signaling and metabolism were downregulated in A2 SMA mESC-derived MNs. Cellular retinol-binding protein (Rbp1) and the cellular RA binding proteins Crabp1 and Crabp2 are expressed in early mouse embryos and may play a role in the development of the CNS [89,90]. The expression of the pluripotency-related protein Zfp42 (Rex1)-whose mRNA was upregulated in A2 SMA mESC-derived MNs, as mentioned earlier-is repressed by RA in F9 teratocarcinoma cells [91]. Rbp1, Crabp1 and Crabp2 mRNA levels were reduced in SMA mESC-derived MNs (Figure 9 and Table S1) as well as in severe SMA spinal cords ( Figure 10). Microarray analysis of PND05 (late-symptomatic) severe SMA spinal cord mRNAs also identify impairments in RA signaling and metabolism [42]. RA regulates many phases during MN differentiation [92]. RA has been implicated in its ability to induce neurogenesis by blocking fibroblast growth factor (FGF) signaling [93]. Furthermore, RA and FGF signaling are sufficient to induce MN differentiation independent of Shh signaling [94]. Table 3. Relationship between changes in mRNA levels measured by RNA-Seq and protein levels measured by 2D-DIGE or immunoblot (*) of selected genes in normal versus SMA mESC-derived motor neurons. The upregulation of pluripotency-related transcripts and the downregulation of transcripts related to neuronal development and activity identified in this study suggest that SMA mESCs may not be differentiating into MNs as efficiently as normal mESCs. The difference between the number of MNs differentiated from A2 SMA and Hb9 control mESCs was not significant (Figure 3). This observation was based on eGFP expression that was driven by the MN promoter HB9. HB9 is a late-stage MN marker [37]. Consistent with the FACS analysis, Hb9 mRNA expression was not significantly altered in SMA mESC-derived MNs even though the mRNA levels for early-stage MN markers like Isl1 and Olig2 were reduced in A2 SMA mESC-derived MNs. The levels of proteins found in MNs-like choline acetyltransferase (ChAT), HB9 and neurofilament-are not altered by SMN deficiency [46]. Taken together, our observations suggest that SMA MNs can initially develop normally but they do exhibit changes in transcripts related to pluripotency and neural development. These transcripts may have novel functions in MNs aside from development.
Microarray analysis of differential gene expression between control and severe SMA spinal cord transcript pools suggest that SMA is a neurodegenerative disease instead of a neurodevelopmental disorder [42]. We did not observe an overrepresentation of apoptosis and cell death transcripts in the pathway and network analyses of our differentially expressed transcriptome data. There is no noticeable cell death in the ventral horn of the spinal cord of severe SMA mouse embryos even though cell death can be detected in the telencephalon [95]. When A2 SMA mESC-derived MNs were plated onto coverslips, we did observe a timedependent loss of cell viability and reduced neurite outgrowth ( Figure 4). Similar reduced neurite outgrowth and cell death are observed in MNs differentiated from induced pluripotent stem cell (iPSC) lines derived from human SMA fibroblasts [96][97][98][99]. Primary MNs cultured from severe SMA mouse embryos exhibit reduced neurite lengths [30]. Knockdown of Smn in zebrafish embryos with morpholino antisense oligonucleotides results in defects in motor axons suggesting early developmental defects [26]. We found that many of the biological pathways downregulated in A2 SMA mESC-derived MNs involved axonal guidance ( Figure 6D). No developmental defects in motor axon formation occur in severe SMA mice [26,95]. In both mouse and fruit fly models for SMA, the maturation and maintenance of neuromuscular junctions is defective and this defect may be the result of denervation injury and/or neurodevelopmental delay [45,[100][101][102][103]. Comparison of presymptomatic control and SMND7 SMA MNs using RNA-Seq reveal deficits in transcripts related to synaptogenesis including agrin (Agrn) [64]. In our isolated mESCderived SMA MNs, we also observed a significant (0.85-fold) decrease in Agrn mRNA levels ( Table S1). Our transcriptome analysis suggests that there may be neurodevelopmental delays in SMA MNs; however, SMA MNs develop normally initially and form connections with target muscles but these connections then atrophy for unknown reasons (reviewed in [21]). Upregulation of pluripotency and cell proliferation transcripts as well downregulation of neuronal development-related transcripts in SMA MNs may be a consequence of denervation and axonal degeneration.
In conclusion, we have identified distinct gene expression patterns in SMA MNs when compared to normal MNs. Pathways upregulated in SMA mESC-derived MNs were involved in pluripotency and cell proliferation whereas common pathways found in the downregulated genes have shown decreases in neuronal markers commonly found in mature and developing neurons. It remains to be determined whether these neuronal marker deficits are a contributing cause or a consequence of the disease. The mechanisms underlying these changes in the transcriptome of SMA MNs will need to be examined in more detail for future studies. Comparison of SMA MN transcriptomes against normal MN RNA transcript profiles will also lead to the identification of novel targets for the development of therapeutics for SMA. Table S1 List of all of the annotated RNA transcripts that are differentially expressed in A2 SMA mESC- Figure 10. Expression of RNA-Seq-identified differentially expressed transcripts in SMA mouse spinal cords. The levels of Crabp1, Crabp2, Isl1, Nkx2.2, Pla2g1b, Smn1 and Vim transcripts were measured in spinal cord total RNAs from control and SMA mice. (A) The magnitude of change (log 2 (fold change)) of selected transcripts in low copy SMN2 SMA (SMN2(89) +/+ ;mSmn 2/2 ; black bars) or high copy SMN2 rescue (SMN2(566) +/+ ;mSmn 2/2 ; grey bars) mouse spinal cord samples (n = 3/genotype) relative to control samples. (B) The magnitude of change of selected transcripts in low copy SMN2 SMA mouse spinal cord samples relative to controls at embryonic day 13.5 (e13.5; black bars) and postnatal day 3 (PND03; grey bars). doi:10.1371/journal.pone.0106818.g010 derived MNs relative to Hb9 control mESC-derived MNs.