Skip to main content

Global transcriptome profile of the developmental principles of in vitro iPSC-to-motor neuron differentiation



Human induced pluripotent stem cells (iPSC) have opened new avenues for regenerative medicine. Consequently, iPSC-derived motor neurons have emerged as potentially viable therapies for spinal cord injuries and neurodegenerative disorders including Amyotrophic Lateral Sclerosis. However, direct clinical application of iPSC bears in itself the risk of tumorigenesis and other unforeseeable genetic or epigenetic abnormalities.


Employing RNA-seq technology, we identified and characterized gene regulatory networks triggered by in vitro chemical reprogramming of iPSC into cells with the molecular features of motor neurons (MNs) whose function in vivo is to innervate effector organs. We present meta-transcriptome signatures of 5 cell types: iPSCs, neural stem cells, motor neuron progenitors, early motor neurons, and mature motor neurons. In strict response to the chemical stimuli, along the MN differentiation axis we observed temporal downregulation of tumor growth factor-β signaling pathway and consistent activation of sonic hedgehog, Wnt/β-catenin, and Notch signaling. Together with gene networks defining neuronal differentiation (neurogenin 2, microtubule-associated protein 2, Pax6, and neuropilin-1), we observed steady accumulation of motor neuron-specific regulatory genes, including Islet-1 and homeobox protein HB9. Interestingly, transcriptome profiling of the differentiation process showed that Ca2+ signaling through cAMP and LPC was downregulated during the conversion of the iPSC to neural stem cells and key regulatory gene activity of the pathway remained inhibited until later stages of motor neuron formation. Pathways shaping the neuronal development and function were well-represented in the early motor neuron cells including, neuroactive ligand-receptor interactions, axon guidance, and the cholinergic synapse formation. A notable hallmark of our in vitro motor neuron maturation in monoculture was the activation of genes encoding G-coupled muscarinic acetylcholine receptors and downregulation of the ionotropic nicotinic acetylcholine receptors expression. We observed the formation of functional neuronal networks as spontaneous oscillations in the extracellular action potentials recorded on multi-electrode array chip after 20 days of differentiation.


Detailed transcriptome profile of each developmental step from iPSC to motor neuron driven by chemical induction provides the guidelines to novel therapeutic approaches in the re-construction efforts of muscle innervation.


Human neuronal tissue lacks regenerative capacity, leaving few treatments available following neuronal injury or neurodegeneration. In the past decade, an interest in direct neuronal reprogramming of stem cells into motor neurons (MNs) has emerged as a solution to generate human neuronal tissue for therapeutic applications. MNs form synapses to potentiate electrical signals from the CNS into peripheral tissues. They play a critical role in the formation of neuromuscular junctions (NMJs), where MN axons terminate on muscle fibers and neurotransmitters are released to trigger muscle contractions. NMJs are cholinergic synapses, where the neurotransmitter acetylcholine (ACh) is released from the presynaptic MN terminal for uptake by postsynaptic ACh receptors on the target muscle cell [1]. This critical function is disrupted in neurodegenerative motor neurons diseases like as Amyotrophic Lateral Sclerosis (ALS).

Several in vitro protocols have been developed to convert progenitor cells, such as human inducible pluripotent stem cells (iPSCs), into MNs [2,3,4]. There are still challenges limiting clinical application of these iPSC-derived MNs. For example, the generation of physiologically active neurons requires a lengthy cell maturation period and often results in a heterogeneous population of neuronal subtypes [5]. Protocol reproducibility can also vary as different cell lineages have unique maturation and functional properties. To address these challenges, we present a 28-day transcriptome study coupled with functional assays. Our main objective was to resolve the underlying mechanisms driving MN differentiation. The results from this study can guide experimental strategies to obtain populations highly enriched with the desired MN subtype.

Here, we followed an established protocol [6] to differentiate iPSCs into MNs using chemically defined media conditions. Efficient neural conversion is based on mimicking in vivo neurogenesis where extrinsic and intrinsic signals are introduced in culture, yielding a relatively pure MN population [7, 8]. During neuronal differentiation, there are cascading effects as signaling pathways activate transcription factors to upregulate expression of MN specific genes. Neural induction of iPSCs is driven by simultaneous inhibition of tumor growth factor-β (TGFβ), activin, Nodal, and bone morphogenic protein (BMP) signaling. Similar to processes that occur during early development, inhibition of those signaling pathways promotes differentiation along the neuronal lineage primarily through inhibition of pluripotency and blocking alternative lineage differentiation. Several other pathways, including the Wnt signaling pathway, regulate neuronal differentiation. The protocol implemented in this study included three core chemical compounds to inhibit TGFβ and BMP signaling pathways and simultaneously activate Wnt signaling. Following neural induction of iPSCs, neuronal progenitor cells were patterned with all-trans retinoic acid (RA), to promote caudal (spinal cord) identity, and ventralization was promoted by activation of Shh signaling with Purmorphamine. Finally, synchronization of the maturation process, through elimination of dividing cells, was aided by inhibition of the Notch signaling pathway resulting in mature MNs.

Genome-wide transcriptome studies provide in-depth knowledge of regulatory pathways that shape cellular morphology and function. Such information is crucial for the design of novel neuron regenerative therapies and cell-based drug discovery platforms. Previous studies based on single cell transcriptomics of Amyotrophic Lateral Sclerosis (ALS) patient-derived iPSCs discovered the underlying mechanisms of disease pathology [9] and the regulatory dynamics of MN differentiation [10]. Others have investigated the iPSC-derived MN axonal transcriptome and found key regulatory pathways presenting potential drug targets for treatment of genetic disorders [11]. A detailed transcriptomics study by Burke et al. proposed an influence of the genetic background of the iPSC donors on each pivotal step of iPSC-initiated corticogenesis [12]. Further, single cell RNA-seq analysis of iPSC-derived spinal MN demonstrated that in vitro differentiation does not produce a homogeneous MN population [8].

Here we performed a comprehensive transcriptomic analysis of in vitro iPSCs-derived MNs to characterize the key principles of MN development by analyzing bulk transcriptome data from five crucial time points of MN neurogenesis: iPSCs (D0), NSCs (D7), post-mitotic MNP cells (D13), eMNs (D18), and mature MNs (D28). We interrogated the transcriptomic signatures using next-generation RNA-sequencing (RNA-Seq) technology and performed in situ validation of key pluripotency and MN specific biomarker expression. Additionally, we characterized the functionality of iPSC-derived MNs via electrophysiological analysis of neuronal network connectivity. Our results corroborate transcriptomic profiles previously identified in the astrocyte-to-neuron transformation process [3]; specifically, upregulation of Shh and Wnt signaling pathways. We observed altered gene expression of TGFβ signaling components during the early stages of iPSC conversion into neural stem cell (NSC) in response to chemical inhibition. However, we detected an upregulation of positive TFGβ regulators in the subsequent steps of neurogenesis.

By applying a LANL-developed Ontology Pathway Analysis software (OPaver) we found that calcium (Ca2+) signaling through cyclic adenosine monophosphate (cAMP) and LPC is downregulated during the conversion of the iPSC to NSC and remains silenced until the final stage of MN maturation when the regulatory pathway becomes a driving force for the neuronal synaptic activity. Another key finding is that during the in vitro development of MN in 2D monocultures, ionotropic nAChR are expressed in the early MN stage and are subsequently replaced by the G protein coupled receptors (GPCR) -type muscarinic AChR in mature MN, probably due to lack of metabolic stimuli that are naturally released by the muscle tissue in vivo. While our findings provide unique insights into the temporal mechanism of iPSC-derived MNs they also indicate the advantages of using a co-culture system, of MN and muscle, to enable the development of physiologically relevant MNs types akin to that obtained in vivo.


A four-step process of iPSC differentiation into MNs

The initiating step for iPSC differentiation to neuronal stem cells (NSC) is driven by inhibition of TGFβ, activin, nodal, and BMP signaling pathways while simultaneously activating the Wnt signaling pathway to sustain cell proliferation. Therefore, to convert undifferentiated iPSCs into a NSC lineage the culture media was supplemented with SB431542 (SB) and DMH-1, inhibitors of activin receptor-like kinases (ALK4, ALK5, ALK7, and BMP), and CHIR99021 (CHIR), a Wnt pathway activator. MNPs and early motor neurons (eMNs) were patterned through activation of the Shh pathway by Purmorphamine (Pur). The final step of MN maturation and specification was aided by a Notch signaling pathway inhibitor, Compound E (CpdE), to synchronize the maturation process through the elimination of dividing cells. Neurotrophic and growth factors were supplemented to facilitate neuronal growth and maturation (Fig. 1a).

Fig. 1
figure 1

Differentiation of WTC-11 induced pluripotent stem cells (iPSC) into motor neurons (MN). a Schematic diagram of the overall experimental design. Shown are the small molecule stimuli for each developmental stage and the corresponding morphological changes in the course of iPSC conversion into MN. The scale bar for the phase contrast images is 400 μm; (b) Representative immunofluorescence images of key pluripotency (Oct4) and pan-neuronal (Nestin, Pax6, and β-III tubulin) biomarkers indicated at efficient iPSC differentiation into neuronal stem cells (NSC) at day 7th post-induction (accumulation of Nestin is shown in orange). Upregulation of Pax6 and β-III tubulin (shown in green) at day 13 of iPSC induction marked the formation of motor neuron progenitor (MNP) cells. c Immunostaining of neuronal structural proteins (beta-III-tubulin in green, MAP2 in orange), motor neuron specific markers (HB9 and ChAT in green), and synaptic vesicles protein (Synaptophysin in green) showed a pure population of mature MN at day 28 post-stimulation of iPSC. d Image quantification of Pax6- and HB9-labeled cells exemplified the percentage of cells that transitioned from NSC to MN. e Immunocytochemistry with anti-ISL1 antibody and quantitative analysis of biomarker positive cells at D28 of the differentiation protocol. The graphs show average of 6 cultures with > 300 cells in random fields for each culture. Scale bar: 100 μm

Within a 13-day period, we observed changes in cell morphology, including a gradual reduction in the cell soma and multiple extensions of thin neurites upon iPSCs conversion into MNPs. Furthermore, we monitored protein accumulation of tissue-specific markers via immunohistochemistry staining (Fig. 1b). Consistent with the morphological features exemplified by high nuclear-to-cytoplasmic ratio and compact multilayer colonies [13], the undifferentiated iPSCs expressed a key pluripotency marker, Oct4 (Fig. 1b; D0). We noted the formation of NSCs at D7 by the disappearance of Oct4 and the accumulation of Nestin and Pax6. Nestin, implicated in cell division and radial axon growth [14, 15], was downregulated by D13 which coincided with increased expression of the pan-neuronal filament protein class III β-Tubulin (βIII-Tub). Cells entered the MNP developmental stage on D13 and transformed to eMNs on D18 when both Nestin and Pax6 expression was significantly downregulated (Fig. 1b, d). Finally, we detected mature MNs featuring high levels of the structural proteins, βIII-Tub and microtubule-associated protein 2 (MAP2); and increased accumulation of MN specific markers: MN homeobox protein (HB9), choline acetyltransferase (ChAT), and the synaptic vesicles protein (Synaptophysin) (Fig. 1c, d). HB9 is an essential transcription factor and early marker of cholinergic neurons [16] and ChAT is an enzyme required for the synthesis of the neurotransmitter acetylcholine (ACh). The changes in morphology were pronounced as dendrites extended and expressed MAP2, a neuron-specific cytoskeletal protein critical for successful projection of dendrites [17, 18]. Accumulation of the pan-MN marker ISL-1 (~ 60% ± 12%, Fig. 1e) together with the HB9 transcription regulator (~ 85% ± 15%, Fig. 1d) at D28 of the iPSC differentiation indicated at formation of cell populations enriched in mature MNs.

MNs form functional network connections

MN activity was characterized by multi-electrode array (MEA) recordings of cellular electrophysiological responses. This non-invasive method allows for repetitive recordings of spontaneous electrical firings at various times during neuronal differentiation. Spontaneous oscillations in extracellular action potential (AP) were recorded after day 20 of neuronal differentiation, and the AP frequency increased the longer the cells differentiated on the MEA (Fig. 2a). By day 31 (D31) the firing patterns became more organized in highly synchronous bursts of network activity. The AP spike rate increased as differentiation progressed (Fig. 2b). The increase in the number of bursts, and percentage of spikes in bursts, signified successful formation of synaptic connections with synchronized AP firings (Fig. 2c).

Fig. 2
figure 2

Electrophysiological analysis of (MNs). a Time-dependent increase of spontaneous action potential (AP) spikes generated by MNs and recorded from a representative electrode (one of 144) on a multi-electrode (MEA) chip. MN were plated on the MEA chip at 18 days post-chemical induction of iPSC (D18) and electrophysiological activity was recorded at the indicated time points of MN maturation (D22-D31) for 1 min. The frequency of AP spikes (b) and the percentage of spikes in the burst (c) dramatically increased within the neuronal networks over time. Shown are the mean and standard error from 12 electrodes at each sampling day post iPSC induction

Unbiased transcriptome analysis of chemically stimulated iPSC differentiation into MNs

We performed whole-genome transcriptome analysis on RNA-Seq platform to explore gene regulatory pathways. Three independent biological replicates were collected and analyzed at each stage: iPSCs, NSCs, MNPs, eMNs, and MNs. The principal component analysis (PCA) of global transcriptome data revealed tight clustering of all replicates, indicating high reproducibility of gene expression profiles for each stage of development. Differentiation timepoints showed distinct genetic programs, as evidenced by the PCA variances (Fig. 3a). Applying pair-wise analyses, we further compared the differential gene expression profiles between cell populations at various developmental stages. The number of differentially expressed genes (DEGs) increased at each timepoint throughout the differentiation process when compared to D0, illustrating the steady transformation of iPSCs into MNs driven by chemical stimuli (Fig. 3b and list of DEG in Additional file 1). The greatest number of DEGs occurred upon the initial culturing of iPSCs in neural differentiation media supplemented with an inhibitor of TGFβ signaling pathway and a Wnt pathway agonist (Fig. 3c; D7-D0, see Additional file 1).

Fig. 3
figure 3

Transcriptome analysis of cellular responses to MN maturation. a Principal component analysis of all samples, n = 3; (b) Histogram of differentially expressed genes (DEGs) (adjusted p < 0.0002, fold change > 3) among D7 - D28 samples in the pairwise comparison with D0; (c) Histogram of DEGs in pairwise comparisons between adjacent timepoints (D0-D7, D7-D13, D13–18, D18–28)

We analyzed RNA-Seq data with the NCBI Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 [19] which categorized DEGs between D0 and D28 by gene ontology (GO) terms (Fig. 4a). Of the upregulated genes, 2118 were recognized by DAVID and 1683 (79.5%) mapped to specific GO terms consistent with a shift from undifferentiated mitotic cells (iPSC) to differentiated cells of the neuronal lineage. The upregulated DEGs were associated with GO terms specific to neuronal development: 98 genes were related to dendrites; 144 genes were involved in cell junctions, 68 genes were associated with axon formation, and 61 genes were involved in the postsynaptic membrane. Of the genes downregulated on D28 versus D0, DAVID recognized 1361 genes of which 1234 (90.7%) mapped to a specific GO. The downregulated genes were characteristic of cells undergoing active DNA replication, transcription, and translation associated with nucleoplasm (440 genes), nucleolus (170 genes), nucleus (556), and cytosol (354) GO terms (Fig. 4a).

Fig. 4
figure 4

Bioinformatics analysis of transcriptomic data. a Gene ontology (GO) analysis of iPSCs and MNs as determined by the NCBI Database for Annotation, Visualization and Integrated Discovery (DAVID). The top four significant GO terms of upregulated and downregulated genes are listed, ranked by p-value, when comparing D0 iPSCs and D28 MNs. b OPaver analysis of 12 KEGG pathways ranked by the number of DEG during the chemical reprogramming of iPSCs into MNs: Pathways in cancer; Neuroactive ligand-receptor interaction; MAPK-signaling pathway; Axon guidance; Ras signaling pathway; PI3K-Akt signaling pathway; Calcium signaling pathway; Rap1 signaling pathway; cAMP signaling pathway; Transcriptional mis-regulation in cancer; Singling pathways regulating pluripotency of stem cells; and Cholinergic synapse

LANL-developed OPaver [20] identified 12 gene signaling pathways that were significantly (p < 0.05) altered during iPSC conversion to MNs (Fig. 4b). The top 5 KEGG pathways represented genes involved in cancer regulation, axon guidance, calcium signaling, PI3K-Akt signaling, and MAPK-signaling pathways (Fig. 4b). The OPaver analysis further highlighted the critical role of membrane proteins in cell-to-cell interactions for neuronal development; specifically, neuroactive ligand-receptor interactions and the development of cholinergic synapses.

By focusing our analysis on cellular development pathways, we found significantly altered gene expression profiles of genes involved in TGFβ, Notch, and Shh signaling pathways (Fig. 5). Contrary to our expectations, 7 days after chemical inhibition of TGFβ, we observed a significant downregulation of genes acting as negative regulators of TGFβ signal transduction: PMEPA1 and TGIF1, repressors of SMAD2 function, and SKIL the nuclear repressor of TGFβ-responsive genes [21]. At the same time, a gene positively regulated by TGFβ signaling, JUNB/AP1, was upregulated (Fig. 5a). These data indicate there was a transient inhibition (lasting less than 7 days) of the TGFβ signaling pathway, when small molecules were present in the culture media to target ALK receptor kinases. In the Notch signaling pathways, we observed a downregulation of Hes1, NOTCH2, and PRKCA after NSCs conversion into MNPs (Fig. 5b; D13). Interestingly, gene expression of Hes5 and JAG1 was upregulated until D18 when cells entered the eMN developmental stage. The final MN maturation phase required the addition of CpdE, a Notch signaling pathway inhibitor [22], along with growth factors: insulin-like growth factor 1 (IGF-1), brain-derived neurotrophic factor (BDNF), and ciliary neurotrophic factor (CNTF). Our results indicated that downregulation of Notch regulators JAG1 and Hes5 gene expression was required for the maturation of MNs (Fig. 5b; D28). Shh signaling is essential for patterning and specification of neuronal cells. Significant activation of a key Shh receptor gene, PTCH1, was observed at D7 prior to the addition of Pur, and remained upregulated throughout the differentiation process but decreased in MNs (D28), indicating that this gene may act as a positive regulator of cell division in human neuronal cells. In contrast, ARHGAP36 and CRMP1 genes were downregulated upon transition of iPSCs to NSCs (D7) and responded to Pur with gradual stimulation of gene expression in MNs at D18–28 (Fig. 5c).

Fig. 5
figure 5

Differential gene expression in key signaling pathways shaping MN differentiation. Gene transcript levels were determined by global RNA-seq analysis of (a) TGFβ, (b) Notch, and (c) Sonic Hedgehog (Shh) pathways and (d) Calcium signaling. Fold change was calculated relative to RNA reads in iPSC (D0). The statistics (average value and standard error) were derived from three independent biological replicas with p < 0.001 determined by R-statistical analysis package

Validation of RNA-Seq data via RT-qPCR

Applying RT-qPCR we further validated a subset of key gene markers of pluripotency, NSCs, MNPs, and MNs which were found to be differentially expressed by RNA-Seq profiling. Both methods detected the rapid decrease of pluripotency marker Oct4 by D7 and increase of Nestin, Pax6, and NgN2 expression reaching peak levels at D13 (Fig. 6a-d). Markers of MN specification (Isl1, HB9, MAP2, and ChAT) increased throughout the course of neuronal maturation (Fig. 6e-f, h). ChAT expression levels gradually increased throughout the differentiation process from NSCs to MNs (Fig. 6h). Collectively, quantitative analysis of gene expression via RNA-Seq closely correlated with RT-qPCR data showing dynamics of gene activities typical for the conversion of pluripotent cells to MNs. We observed transient activation of Nestin, a key regulator of cytoskeletal dynamics and cell division in NSC, and stable expression of pan-neuronal markers Pax6, NgN2, and MAP2. Our data indicated there was gradual accumulation of gene transcripts responsible for the synthesis of neurotransmitter, ChAT, and transcription factors, Isl1 and HB9, whose concerted actions shape the phenotype and physiological activity of MNs.

Fig. 6
figure 6

Comparative analysis of differential gene expression of key tissue development markers by RT-qPCR (gray) and RNA-seq (black). Pluripotency and NSC markers validated include (a) Oct4, transcription factor that maintains self-renewal and pluripotency; (b) Nestin, a filament protein marker of neural stem cells; (c) Pax6, a transcription factor that drives neurogenesis; and (d) NgN2, a neuronal-specific transcription factor. Motor neuron specification markers validated include (e) Isl1, a transcription factor required for motor neuron generation; (f) Map2, a neuron-specific cytoskeletal protein; (g) HB9, an early marker of cholinergic neurons; and (h) ChAT, an enzyme required for acetylcholine synthesis. Shown are the averages and standard error from three independent biological replicas from the iPSC to MN differentiation trajectory. The RNA-Seq transcripts were normalized to the total read per analyzed sample (in FPKM: fragments per kilobase per million mapped fragments) and the transcript levels determined by RT-qPCR were normalized to GAPDH as the endogenous sample control. Fold change was calculated for each developmental stage relative to transcript levels in iPSC (D0). Statistical significance (p < 1.5e-6) was determined with Student t-test

Cell signaling pathways driving iPSC conversion to MNs

The GO analysis, conducted by OPaver, identified 12 pathways that included a high number of DEGs as iPSCs were differentiated into MNs. We further investigated the DEG pattern of individual genes in the pathways listed below which had a critical role in shaping neuronal development and function (Tables 1, 2, 3, 4, 5). In each description below, ‘upregulated’ and ‘downregulated’ gene expression refers to a significance of p < 0.05 and log2 fold change ≥2 between sampling timepoints.

Table 1 DEGs identified in KEGG pathway: MAPK-signaling. Differential expression indicates the average log2-fold change in RNA-Seq transcript levels from three independent experiments at each sampling period: D7 vs D0, D13 vs D7, D18 v D13, and D28 v D18. The p-values were ≤ 0.001. Negative log2-fold change corresponds to gene downregulation and the positive values indicate gene activation
Table 2 DEGs identified in KEGG pathway: Calcium signaling. Differential expression indicates the average log2-fold change in RNA-Seq transcript levels from three independent experiments at each sampling period: D7 vs D0, D13 vs D7, D18 v D13, and D28 v D18. The p-values were ≤ 0.001. Negative log2-fold change corresponds to gene downregulation and the positive values indicate gene activation
Table 3 DEGs identified in KEGG pathway: Neuroactive ligand-receptor interaction. Differential expression indicates the average log2-fold change in RNA-Seq transcript levels from three independent experiments at each sampling period: D7 vs D0, D13 vs D7, D18 v D13, and D28 v D18. The p-values were ≤ 0.001. Negative log2-fold change corresponds to gene downregulation and the positive values indicate gene activation
Table 4 DEGs identified in KEGG pathway: Axon guidance. Differential expression indicates the average log2-fold change in RNA-Seq transcript levels from three independent experiments at each sampling period: D7 vs D0, D13 vs D7, D18 v D13, and D28 v D18. The p-values were ≤ 0.001. Negative log2-fold change corresponds to gene downregulation and the positive values indicate gene activation
Table 5 DEGs identified in KEGG pathway: Cholinergic synapse. Differential expression indicates the average log2-fold change in RNA-Seq transcript levels from three independent experiments at each sampling period: D7 vs D0, D13 vs D7, D18 v D13, and D28 v D18. The p-values were ≤ 0.001. Negative log2-fold change corresponds to gene downregulation and the positive values indicate gene activation

MAPK-signaling pathway

MAPK-signaling (Table 1) has an important role in mediating neuronal differentiation and survival [23]. Consistent with our expectations, the transition from iPSCs to NSCs was guided by inhibition of signal transduction pathways supporting self-renewal and activation of genes regulating neuron-specific development pathways. Gene expression of Neurotrophic Tyrosine Kinase Receptor type 2, NTRK2, was prominently upregulated at D7 when the NSC population was established, while a subset of genes encoding for γ-subunits of voltage-dependent calcium-channels (Cavs) were dramatically downregulated (CACNG5, CACNG7, and CACNG8), along with MYC, Ras guanyl-releasing protein 2 (RASGRP2), and GNG12. Gene groups encoding the auxiliary subunits of high-voltage activated (HVA) Cavs were upregulated during the early stages of neuronal development [24]: CACNB3 (induced 4-fold at D7) and CACNA2D3 (16-fold activation at D7-D13). The late stage of MN formation (D18-D28) was marked by upregulation of CACNA1E and CACNG2. Several gene groups encoding auxiliary subunits of low-voltage activated (LVA) Cavs (CACNG5 and CACNG8), were initially downregulated in NSC and later upregulated during the MN maturation stage (D18-D28). Accumulation of Cavs gene transcripts throughout the MN differentiation process is consistent with the increased excitability of neuronal cells (Fig. 2).

Other key components of MAPK-signaling showed oscillatory patterns of gene expression throughout MN formation and maturation. Gene transcripts of KIT and its ligand, KITLG, which function in stem cell maintenance, were upregulated in NSCs (D7) and downregulated during the transition from MNP to eMN (D13-D18). KITLG transcription was again upregulated in MNs (D18–28). KIT signaling promotes cell survival through the activation of PIK3, PLC, and AKT1 pathways and we detected upregulation of KITLG gene expression in MNs which may be a response to the addition of IGF-1 in the culture media. Upregulation of various genes regulating cell survival and proliferation in the post-mitotic MN, such as KDR, PTPRR, PTPN5, and PRKCB, further exemplifies the MN response to IGF-1 pro-survival stimulus.

Calcium (Ca2+) signaling pathway

Ca2+-signaling (Table 2) is a critical component of synaptic activity and neuronal function [25]. Our RNA-Seq data indicated the conversion of iPSCs into NSCs is marked by a significant downregulation of a gene encoding for a serotonin-specific GPCR, 5-hydroxytryptamine receptor 7 (HTR7). In addition, adenylate cyclase 2 (ADCY2) was also downregulated at this time and is known to act immediately downstream of HTRs. Together, the downregulation of HTR7 and ADCY2 indicates a functional silencing of the cAMP-dependent signaling pathway. Furthermore, multiple genes involved in the activation of Ca2+ signaling were downregulated: GPCR subunit alpha-14 (GNA14), an activator of phospholipase C (LPC); two gene groups encoding inositol 1,4,5-triphosphate receptors (ITPR2 and ITPR3); and ryanodine receptor 2 (RYR2). (Fig. 5d and Additional File 2) One GPCR from the chemokine family, CXCR4, involved in the cytosolic Ca2+ mobilization and cell migration, was also markedly upregulated (16-fold) during the transition into NSCs.

A hallmark of the NSC to MNP transition was the accumulation of receptor tyrosine-protein kinase erB-4 (ERBB4) transcripts. ERBB4 is activated by epidermal growth factor proteins, neuregulins2/3, to shape the development of neuronal cells through activation of MAPK3/ERK [26, 27]. ERBB4 expression was upregulated by 18-fold in D7-D13, followed by upregulation of CACNA1B, CACNA1E, GRIN2A and P2RX3 (D13-D18). Formation of eMNs (D13-D18) coincided with upregulated gene expression of HVA N- and R-type (CACNA1B and CACNA1E) voltage-dependent Ca2+ channels, glutamate ionotropic receptor (GRIN2A), and the ligand-gated ion channel responsible for peripheral pain responses, P2RX3. These gene expression profiles indicate the selective pressure of chemical stimuli lead to the formation of eMNs with characteristics of the sympathetic nerve system [28].

In the final stage of MN maturation (D18-D28), we observed upregulation of gene groups activating Ca2+ signaling through the cAMP-dependent pathway: HTR7 and ADCY2. We also detected upregulated expression of the intracellular Ca2+ ryanodine receptor 2 (RYR2), signifying cell readiness to release Ca2+ from the sarcoplasmic reticulum in response to adrenergic (sympathetic) stimulation (Fig. 5d). This finding was consistent with the upregulation of ADRA1A, which encodes the adrenergic receptor alpha subunit 1α. Furthermore, gene activation of both ionotropic (GRIN1 and GRIN2D) and metabotropic (GRM1) glutamate receptors, together with the upregulation of neuromedin-K receptor (TACR3), accentuated the role of the IP3-Ca2+ second messenger system to activate ERK1/2 as well as classical PKC signaling pathways [29,30,31]. We also detected higher levels of PRKCB gene transcripts in MNs compared to eMNs (Table 2).

Interestingly, G-coupled muscarinic acetylcholine receptors (mAChRs) were upregulated in mature MNs (D28). Genes encoding for excitatory (CHRM3) and inhibitory (CHRM2) mAChRs were specifically activated in late stage MN maturation (Table 2, D18-D28). The excitatory mAChRs stimulate PLC, triggering IP3 and diacylglycerol signaling pathways (Additional file 2). Since mAChRs are present on parasympathetic neurons [32], our data suggest that the iPSC-induced MNs bear the excitatory characteristics of parasympathetic neurons.

Neuroactive ligand-receptor interactions

The majority of neuronal receptors are GPCRs and they function to regulate signal transduction pathways and shape cellular physiological responses. Our genome wide transcriptome analysis revealed that iPSC-induced MNs differentiation is marked by the expression of neuroactive receptors which are responsive to glutamate, ACh, and catecholamines, orexin, and prostaglandins (Table 3).

During the first step of iPSC commitment to a neuronal lineage, multiple glutamate receptors were altered, notably, gene expression levels of two ionotropic glutamate receptors were upregulated in the NSCs, GRIA1 and GRIK1, while their paralog genes, GRIA4 and GRIK4, were downregulated. This suggests a specific requirement for type 1 ionotropic glutamate receptors early in neuronal development. Neuromedin and apelin receptors (NMUR2 and APLNR) were upregulated in NSCs, as well as GPCR response to sphingosine-1-phosphate signaling receptor (S1PR1), suggesting changes in cytoskeleton dynamics and mitosis [33]. Genes encoding for nicotinic acetylcholine receptors (nAChR) and cannabinoid receptors (CNR1and CHRNB4) were upregulated during the conversion of NSCs into MNPs and remained activated until eMNs were formed.

Formation of eMNs was marked by reactivation of genes encoding glutamate receptors which were downregulated at the NSC stage. Gene activity of various ionotropic glutamate receptors peaked when eMNs became mature MNs (D18-D28). Activation of genes regulating the synthesis of various neurotransmitters defined the transition from eMNs to MNs. The transcript levels of the neuropeptide tachykinin, and its receptor (TACR3), were upregulated in MNs. In spinal neurons, tachykinins evoke synthesis and release of ACh, histamine, catecholamines, and GABA. We also detected differential expression of the somatostatin receptor (SSTR1) indicating that the mature MN population can secrete neurotransmitters [34]. The glucocorticoid receptor (NR3C1) was upregulated in MNPs and MNs and two new classes of neuroactive receptors were upregulated in the last stage of MN maturation: hypocretin/orexin (HCRTR2) and prostaglandin E receptor 4 (PTGER4). Our transcriptomic data indicate that the replacement of ionotropic nAChR (CHRNA3) and dopamine receptor (DRD2) expression with adrenergic receptor (ADRA1A) and mAChRs (CHRM2 and CHRM3) marked the formation of mature MNs. The controversial activation of mAChRs and genes encoding excitatory neurotransmitters, including glutamate and catecholamines (ADRA1A), could stem from the formation of a heterogeneous population consisting of parasympathetic and sympathetic MNs.

Axon guidance and cholinergic synapse

Pathways regulating axon guidance and synapse development were upregulated during iPSC-to-MN differentiation (Tables 4 and 5, respectively). These functions are critical for neuronal development, maintenance, and repair mechanisms.

We found netrin-G1 (NTNG1) to be initially downregulated as iPSCs transitioned into NSC, but then upregulated in eMNs and MNs from D13–28. Deleted in colorectal carcinoma (DCC) was also significantly upregulated as iPSCs transition through NSCs to become MNPs, indicating its critical role in prompting axon extension. SLIT1 was upregulated at D7 and both SLIT1 and SLIT2 were upregulated as the NSCs became MNPs (D7-D13). HB9 tightly regulates Robo2 expression to regulate motor axon guidance in ventrally-projecting MNs [35] and we observed ROBO2 downregulation in NSCs, during which time HB9 transcripts were initially accumulating (Fig. 6g).

In vivo as MN axons extend, they have the potential to terminate on a muscle fiber at a cholinergic synapse, the site of ACh neurotransmission. In the pre-synaptic neuron, ACh synthesis is driven by choline O-acetyltransferase (ChAT) and we observed ChAT gene expression to be increased in both eMNs and MNs. Expression of ion-channel-coupled nAChR subunits were also significantly altered: CHRNA3 was regulated in NSCs but downregulated in MNs and CHRNA4 was upregulated at the eMN stage. In the postsynaptic membrane these receptor components would contribute to neuroactive ligand-receptor interactions as previously described. Gene expression of acetylcholinesterase (AChE) was significantly increased late in MNs maturation (D18-D28), which encodes an enzyme responsible for the hydrolysis of ACh into choline and acetic acid, a critical step in ACh recycling in the cholinergic synapse. Membrane transporter solute carrier family 5 member 7 (SLC5A7) is then able to import choline back into a cholinergic neuron for subsequent ACh synthesis. We found gene expression of SLC5A7 to be increased in both the eMNs and mature MNs.

Heterogeneity of MN populations

As previously reported [8], in vitro differentiation of human iPSC resulted in heterogeneous populations consisting of several MN subtypes. Applying the same differentiation protocol [6] and quantitative image analysis (Fig. 1d-e), we observed that two pan-MN transcription regulators ISL-1and HB9 were present in ~ 60 (±12%) and ~ 85 (±15%) of the cell populations, respectively (Fig. 1d-e). Thus, our quantitative image analysis data is in agreement with the 68 ± 7% double positive ISL1+/HB9+ MN populations reported by Thiry et al [8]. Based on our image analysis on single cell level (Fig. 1) and quantitative RT-PCR of bulk transcriptome (Fig. 6) we could infer that iPSC-derived neuronal populations are highly enriched on MN in general although without detailed analysis of the molecular markers characteristic to each MN subtype we are not able to conclude which one is predominant. Combination of transcription master regulators' activity determines the MN functional specifications (Table 6). Based on the quantitative image analysis, majority of the MNs were positive for the HB9 gene (Fig. 1c-d), which remains activated in mature MNs from the median and hypaxial motor columns innervating the long muscles of the back and the respiratory muscles, respectively, while HB9 is downregulated in the medial neurons of the lateral motor column (LMC) innervating the limb muscles. The meta-transcriptome data from our study showed significant (p < 1e-5) activation (log2 > 2) of key molecular markers specific to various MN subtypes ranging from the preganglionic column (SMAD1 and ZEB2), spinal accessory column (RUNX1 and PHOX2B), and the lateral LMC (FOXP2 and ALDH1A2). Taken together, our results (summarized in Table 6) strongly indicate that in vitro derived MN populations are heterogeneous consisting of several functional subtypes.

Table 6 MN specification based on log2-fold change in RNA-Seq transcript levels at D28 versus D0 of the differentiation protocol. The p-values were calculated with the R-package. Negative log2-fold change corresponds to gene downregulation and the positive values indicate gene activation. indicates the gene is expected to be downregulated in the MN subtype


This comprehensive study documents transcriptomic and morphological changes of iPSCs as they differentiate into motor neurons (MNs) in vitro. We performed RNA-Seq meta-transcriptome analysis of human iPSC and the four cell types representing stages of spinal motor neuron differentiation (NSC, MNP, and eMN) and maturation (MN). The results from our genome-wide transcriptome study revealed basic developmental principles of in vitro neurogenesis from iPSC that have not been elucidated by previous studies while confirming the regulatory role of TGFβ, Notch, and Shh signaling pathways in the formation of adult spinal motor neurons. We further corroborated the findings from the next-generation RNA-Seq analysis with RT-qPCR gene expression assays and immunohistochemistry profiles of key pluripotency and MN markers. Applying novel OPaver software we found strict temporal correlation between the formation of functional neuronal network connections on a MEA chip and the expression of genes involved in the regulation of Ca2+ signaling. We found that cAMP-regulated Ca2+ signaling was inhibited on gene transcript level when iPSC traversed through the neurodevelopmental stages and was reactivated in the final stage (D28) of MN maturation when we were able to detect neuronal synaptic activity via recordings of spontaneous AP firings.

Inhibition of Ca2+ signaling is required for the transition from iPSC to NSC

Intracellular Ca2+ levels, known as calcium transients, regulate neuronal subtype and neurotransmitter specifications [36, 37]. Our global transcriptome data showed activation of genes regulating Ca2+ entry, specifically voltage-sensitive calcium channels and ionotropic glutamate receptors, not sooner than at the eMN developmental stage (D18 vs D13, Table 1 and Fig. 5d). Ca2+ signaling comprises a cascade of molecular interactions and biophysical events, which translate extracellular signals to intracellular responses via increase of cytoplasmic Ca2+. This can be activated by neurotransmitters, hormones, and growth factors, chemical and electrical stimuli, causing membrane excitation. Two fundamental mechanisms regulate Ca2+entry through protein channels: voltage-dependent Ca2+ channels and ligand-gated channels. The latter are highly diverse, non-Ca2+-specific, and greatly represented by the family of guanine-binding GPCR [38]. Significant upregulation in differentiated MN of high-voltage activated Cavs genes, such as CACNB3, is responsible for tight control on intracellular Ca2+ signaling through regulation of calcium (Ca2+) entry and direct interaction with phospholipase C-coupled (PLC) and inositol trisphosphate (IP3) receptors [39]. Upregulation of CXCR4 by 16-fold at D0-D7 underlines the importance of intracellular Ca2+ signaling inhibition as the neuronal cell lineage commitment process begins. The CXCR4/CXCR7 complex recruits β-arrestin to trigger the canonical GPCR pathway activating ERK1/2, p38, and SAPK, while inhibiting both Ca2+ mobilization and cAMP signaling [40]. Upregulation of neuromedin and apelin receptors (NMUR2 and APLNR) in NSCs suggest a role for phosphoinositide (PI) signaling pathways that inhibit cAMP production upon Ca2+ mobilization and possibly in regulating cytoskeleton dynamics, cell growth, and hormone release [41, 42].

Our findings underline the role of temporal gene regulation of Ca2+ signaling in motor neuron development in vitro. It has been previously demonstrated that voltage-gated Ca2+ influx activated by cAMP is instrumental in the maturation of neuronal progenitor cells into functional neurons [43]. The comparative transcriptome study presented here revealed that genes regulating cAMP synthesis (ADCY2), voltage-gated calcium channels (CACNA1B and CACNG5), and the receptor regulating the intracellular Ca2+ homeostasis (RYR2) were immediately inhibited upon transition of iPSC into NSC and were re-activated in mature MN. Collectively, our data suggests that inhibition of pathways regulating the Ca2+ transients is required for the successful transition from pluripotency to neuronal progenitor cells. Our findings corroborate previous studies that have demonstrated that low cellular excitability is vital for cell migration while increase in the Ca2+ transients stops migration and promotes dendrite formation in cortical neurons [44].

MN specification in vitro is driven by the chemical stimuli in growth media and bare the characteristics of the parasympathetic nervous system

The transcriptome analysis revealed the underlying molecular changes in differentiation of iPSC to mature MNs. Specifically, we observed that short-term inhibition of TGFβ signaling was sufficient to push iPSCs into NSCs and that continuous activation of Notch and Shh signaling pathways ensured morphogenesis and cell survival throughout the MN differentiation and maturation process. While the majority of changes were characteristic of neuronal development, deeper analysis revealed the involvement of more unexpected genes as well as unique temporal changes. For example, the oscillatory pattern of SKIL gene expression throughout the neuronal cell differentiation process suggests a regulatory feed-back loop that balances survival and cell differentiation programs. PTCH1, a Shh receptor gene, is consistently and significantly upregulated from D7 to D28 suggesting a role as a positive regulator in neuronal cell division. Binding of Shh to PTCH1 results in the release of the smoothened protein initiating cell proliferation. Upon ligand binding, PTCH1 is trafficked away from the Shh positive regulator, G-coupled receptor SMO, resulting in downstream signal transduction [45]. Downregulation of PTCH1 by D28 supports the possible inhibition of cell proliferation and subsequent differentiation to neurons. ARHGAP36 and CRMP1, both of which are downregulated simultaneously, encode for Rho GTPase activating protein 36 and collapsin response mediator protein 1, are implicated in semaphorin-induced growth cone and axon guidance [46, 47] and may regulate cell division and morphogenesis of MNs.

Applying quantitative gene transcript analysis and immunohistochemistry, we observed that ChAT, a marker of mature neurons, was upregulated as soon as iPSC committed to a neuronal lineage (D7). In addition, the well-orchestrated expression of key genes such as Nestin, Pax6, NgN2, MAP2, Isl1 and HB9 corroborated the phenotypic changes to MNs. The patterns of gene expression reported here are consistent with similar studies using small molecules to drive neuron reprogramming from astrocytes [3] or fibroblasts [48].

While Shh drives motor neuron formation over intermediate neurons, there are multiple subtypes of MNs. In the final stage of MN maturation, we observed that genes encoding ionotropic nAChR (CHRNA3) and dopamine receptor (DRD2) were downregulated while adrenergic receptor (ADRA1A) and mAChRs (CHRM2 and CHRM3) gene expression was activated. Since mAChRs are found exclusively on the neurons from the parasympathetic system [49] our data indicates that in vitro iPSC-derived MN in monoculture have the characteristics of the parasympathetic neurons of the peripheral nervous system.

In addition to the parasympathetic neuron-specific mAChRs gene activation, genes encoding excitatory neurotransmitters specific to the sympathetic neurons, including glutamate and catecholamines (ADRA1A), were significantly upregulated in the population of mature motor neurons. This could stem from the formation of a heterogeneous population consisting of parasympathetic and sympathetic MNs. Previous studies based on single-cell RNA-Seq analysis of iPSC-derived spinal MN have demonstrated that the protocol for in vitro differentiation produces a mixed population of MN subtypes with a predominant (58%) sub-population of lateral MNs and several minority sub-populations, including hypaxial motor column (19%) and median motor column (6%) MNs [8].

The neuroactive ligands define electrophysiological activity of in vitro iPSC-derived mature MNs

The electrophysiological activity of mature MN reflects their ability to form functional neuronal network connections. Neuronal precursors and neurons are capable of spontaneous electrical activity [50], and we observed increased electrical activity as synapses and neural circuits formed. As iPSCs differentiated after plating on the MEAs, dendrite projections formed as axons stretched between neuronal cell bodies. These networks of MNs created numerous synapses to propagate nerve impulses. The MEAs were able to capture spontaneous firings of extracellular action potentials as single spikes of activity early in differentiation which progressed to more frequent firing bursts as the iPSC-derived MNs matured, demonstrating peak activity on D31. While the length and amplitude of an action potential are always the same, an increase in the stimulus caused an increase in the frequency of an action potential indicative of an enhanced response.

MNs are often characterized by their role in forming cholinergic synapses at NMJs. MN axons terminate on muscle fibers and nerve impulses are translated into muscle contractions as the neurotransmitter ACh is released from the presynaptic MN terminal for uptake by postsynaptic ACh receptors on the target muscle cell [1]. Understanding the underlying mechanism of cholinergic development provides insights about potential treatments for neurodegenerative diseases and strategies to develop countermeasures to chemical toxicology exposure. AChE inhibitors are used to treat Alzheimer’s and Parkinson’s diseases [51, 52] and cholinergic agonist treatments have shown to improve memory function [53]. Similar strategies have been employed to treat organophosphate-induced neurotoxicity [54, 55]. Here we show how culture conditions imprint distinct fates and future efforts to co-culture of MNs with other cell types, such as muscle, may cause different neuronal specification. Surprisingly, we observed the expression of nAChRs typically found on the NMJ to be downregulated in mature MN monocultures while the mAChRs were upregulated. These results can further be applied to characterize aberrant neuronal functions following neurodegeneration or exposure to chemical toxins.


This study revealed significant changes to the transcriptome, protein expression, and electrical function as iPSCs differentiated into mature MNs. To ensure cellular mobility that is essential for tissue layer formation, genes regulating Ca2+ signaling are downregulated throughout the cell differentiation process and are activated in mature MN. Differentiation of iPSC into MN in vitro monocultures leads to formation of both parasympathetic and sympathetic neuronal types. Downregulation of tumor suppressing genes upon chemical conversion of iPSC to mature MN underlines the danger from direct application of the technology in vivo. Understanding the underlying molecular and cellular cues involved in MN differentiation of iPSCs has the potential to enable the discovery of novel treatments for neural injuries.


Culturing and differentiating iPSCs

Human iPSCs (WTC-11, Coriell Institute) were cultured and maintained on vitronectin (ThermoFisher Scientific) treated culture plates in Essential8 medium (ThermoFisher Scientific). Differentiation of iPSCs into MNs was directed as previously described [6] with slight modifications that substituted the NSC growth in suspension (D7-D13) with a growth on vitronectin substrate to improve the cell survival rate and to ensure gene expression profile characteristic to motor neuron progenitor cells (MNP). Briefly, iPSCs were cultured in neural media which consisted of 1:1 DMEM/F12 and Neurobasal medium supplemented with N2, B27, 1x Glutamax and 1x penicillin/streptomycin (all from ThermoFisher Scientific), and 0.1 mM ascorbic acid (StemCell Technology). On days 0–6 3 μM CHIR99021 (StemCell Technology), 2 μM DMH-1 (Tocris) and 2 μM SB431542 (StemCell Technology) were added to the neural medium; days 6–12 the same media was used with the addition of 0.1 μM RA (StemCell Technology) and 0.5 μM Pur (Sigma); days 12–18 cells were maintained with 0.1 μM RA and 0.5 μM Pur added to the neural media; finally from day 18 on cells were cultured with 0.5 μM RA, 0.1 μM Pur and 0.1 μM CpdE (StemCell Technology), IGF-1, BDNF, and CNTF (all from R&D Systems, 10 ng/ml each). For optimum neuronal network formation, cells at the MNP differentiation stage (D12) were plated on MEA at 60% area density.

RNA extraction

Cells were lifted with accutase, pelleted by centrifugation, and stored at − 20 °C. Total RNA was extracted from cell pellets using RNeasy Mini Kit (Qiagen), following the recommendations of the manufacturer. After DNase digestion by Turbo DNA-free kit (ThermoFisher Scientific), samples were quantified and divided for qPCR and transcriptomic analyses.

RNA sequencing analysis

Extracted and DNase-treated RNA was quantified using the Qubit 4 Fluorometer (ThermoFisher) with the High Sensitivity RNA reagents and Bioanalyzer (Agilent) with RNA 6000 Pico reagents. Ribosomal depletion, DNA conversion, and library preparation was performed on all samples using the Illumina TruSeq Stranded Total RNA kit. 151 base pair reads were sequenced on the Illumina NextSeq. Across fifteen samples (three independent experiments x five time points) the total number of reads generated for each sample ranged from approximately 26 million to 40 million reads. Sequencing data was quality trimmed using FaQCs [56] with a quality score cutoff of Q20. Differential expression analysis was performed using PiReT [57] V 0.3.2 and utilizing DEseq2 [58] default parameters and setting a q-value of 0.05 (false discovery rate metric). The experimental design file (provided in the supplementary material) was used to dictate the replicate sample ID’s and sequencing data to be used in the PiReT analysis. Human genome version hg38 was used as the reference genome. KEGG [59, 60] pathway mapping was performed using Omics Pathway Viewer - ‘OPaver’ (Li, unpublished). Raw RNA-Seq reads were deposited in the NCBI SRA database under the accession numbers SRR11994167- SRR11994181. Metadata for each sample are also accessible under NCBI BioProject PRJNA638768.

Quantitative reverse transcription PCR (RT-qPCR)

Three independent experiments were run in duplicate using a 7500 Fast Real-Time PCR System (Applied Biosciences). Fifty ng of each RNA sample were probed for motor neuron differentiation markers using Taqman RNA-to-CT 1-Step Kit (Applied Bioscience) in a 25 μl volume according to the manufacturer’s instructions. Taqman probes included NEUROG2 (Hs00935087_g1), ChAT (Hs00758143_m1), Isl1 (Hs01099686_m1), PAX6 (Hs01088114_m1), MAP2 (Hs00258900_m1), Nestin (Hs04187831_g1), Oct4 (Hs00999632_g1), and HB9 (Hs00907365_m1). Two endogenous controls, actin (Hs99999903_m1) and GAPDH (Hs01922876_u1), were analyzed by RT-qPCR and no significant differences were observed in their expression levels (data not shown); thus, all data shown are normalized to GAPDH.

Immunocytochemistry staining and analysis

Nunc Lab-Tek chamber slides (ThermoFisher Scientific) were coated with vitronectin, seeded, and subsequently fixed with 4% paraformaldehyde, permeabilized with 0.4% Triton X-100, blocked with 3% BSA in PBS for at least 1 h. Samples were incubated overnight at 4 °C with primary antibody solutions (Table 6) diluted in PBS containing Image-iT FX signal enhancer. Cells were washed with PBS three times prior to incubation with NucBlue Fixed Cell Reagent, Image-iT FX signal enhancer, and secondary Alexa-488-, Alexa 555-, or Alexa 647-conjugated antibodies at 37 °C for 2 h (1:1000 dilution, ThermoFisher Scientific and Jackson ImmunoResearch). Table 7 summarizes the antibodies and their concentrations applied in this study. After three PBS washes, the media chambers were removed from the glass slide, coverslips were mounted using ProLong Diamond Antifade Mountant and cells were examined using fluorescence microscopy (Zeiss Observer Z.1). For biomarker quantification, images were acquired with an AxioCam camera connected to Axio Observer Z1 microscope, using ZEN software. For each culture, images were taken with a 20X objective choosing fields with > 100 cells. Images from 5 random fields per culture condition were analyzed with Image J with the following parameters: (i) manual threshold, (ii) fill holes, (iii) enhance contrast and saturated pixels at 0.3%, (iv) watershed separation, and (v) particle analysis at size 100-infinity and circularity at 0.2–0.1. The fraction of cells expressing biomarker proteins was calculated as percent of total cells labeled with DAPI.

Table 7 List of primary antibodies

Functional analysis of MNs on microelectrode array (MEA)

Cells were seeded on MEA chips, either in 60MEA200/30iR-Ti arrays or 24-well Plate with PEDOT Electrodes on Glass, 24 W300/30G-288 (Multichannel Systems). MEA’s were coated with poly-D-lysine and vitronectin. Before recording, the MEA chips were moved to the MEA2100 system (MultiChannel Systems) equipped with temperature control and allowed to equilibrate for 10 min before recording. The data were acquired using Multi Channel Experimenter or Multiwell Screen (MultiChannel Systems) at a sampling rate of 20 kHz for 2 min at 37 °C. Data were filtered using Butterworth band pass filter with 200 Hz cutoff frequency and threshold of 5 x SD were set to minimize both false-positive and missed detection. The representative electrodes were selected for analysis of mean spike frequency and percentage of spikes in the burst.

Transcriptome analysis

DAVID Functional Annotation Bioinformatics Microarray Analysis can be accessed at The raw transcriptomic data of D0-D28 significant DEGs (p < 0.05) included 2242 upregulated and 1438 downregulated terms, available in the supplemental material. Each list of Homo sapiens genes was independently analyzed by DAVID, to generate an analysis of associated gene ontology (GO) terms.

OPaver (Li, unpublished), is a web-based tool to integrate multiple types (e.g. transcriptomics, proteomics and metabolomics) and time series of data to KEGG biochemical pathways maps. This software analysis tool was developed at Los Alamos National Laboratory. In this case, OPaver was utilized to map significantly differentially expressed genes (p < 0.05) identified in the DEseq2 analysis performed by PiRet. Differential expression calculation from DEseq2 in Log2 fold change and associated genes (provided in the supplementary material) were used as input for the OPaver software.

Availability of data and materials

The datasets generated for this study can be found in the NCBI SRA database under accession numbers SRR11994167- SRR11994181, NBCI BioProject PRJNA638768, [] and NCBI BioSamples SAMN15207814-SAMN15207828.





Amyotrophic lateral sclerosis


Action potentials


Brain-derived neurotrophic factor


Class III β-tubulin


Bone morphogenic protein


cyclic adenosine monophosphate


Voltage-dependent calcium-channels


Choline acetyltransferase




Ciliary neurotrophic factor


Compound E


Day 0


Day 7


Day 13


Day 18


Day 28


Day 31


Database for annotation, visualization and integrated discovery


Differentially expressed genes


Early motor neurons


Gene ontology


G protein-coupled receptors


High-voltage activated


Inositol trisphosphate


Induced pluripotent stem cells


Insulin-like growth factor 1




Low-voltage activated


Los Alamos national laboratory


Muscarinic acetylcholine receptors

MAP 2:

Microtubule-associated protein 2


Multi-electrode array


Motor neurons


Motor neuron progenitors


Nicotinic acetylcholine receptor


Neuromuscular junctions


Neural stem cells


Ontology pathway analysis software


Principal component analysis


Phospholipase C




All-trans retinoic acid






Sonic hedgehog.


Tumor growth factor-β


  1. Legay C, Mei L. Moving forward with the neuromuscular junction. J Neurochem. 2017;142(Suppl 2):59–63.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  2. Sances S, Bruijn LI, Chandran S, Eggan K, Ho R, Klim JR, et al. Modeling ALS with motor neurons derived from human induced pluripotent stem cells. Nat Neurosci. 2016;19(4):542–53.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Ma NX, Yin JC, Chen G. Transcriptome analysis of small molecule-mediated astrocyte-to-neuron reprogramming. Front Cell Dev Biol. 2019;7:82.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Davis-Dusenbery BN, Williams LA, Klim JR, Eggan K. How to make spinal motor neurons. Development. 2014;141(3):491–501.

    Article  CAS  PubMed  Google Scholar 

  5. Amoroso MW, Croft GF, Williams DJ, O'Keeffe S, Carrasco MA, Davis AR, et al. Accelerated high-yield generation of limb-innervating motor neurons from human stem cells. J Neurosci. 2013;33(2):574–86.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Du ZW, Chen H, Liu H, Lu J, Qian K, Huang CL, et al. Generation and expansion of highly pure motor neuron progenitors from human pluripotent stem cells. Nat Commun. 2015;6:6626.

    Article  CAS  PubMed  Google Scholar 

  7. Petros TJ, Tyson JA, Anderson SA. Pluripotent stem cells for the study of CNS development. Front Mol Neurosci. 2011;4:30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Thiry L, Hamel R, Pluchino S, Durcan T, Stifani S. Characterization of Human iPSC-derived Spinal Motor Neurons by Single-cell RNA Sequencing. Neuroscience. 2020;450:57–70.

  9. Namboori SC, Thomas P, Ames R, Garrett LO, Willis CR, Stanton LW, Bhinge A. Single cell transcriptomics identifies master regulators of dysfunctional pathways in SOD1 ALS motor neurons. bioRxiv. 2019:593129.

  10. Sagner A, Gaber ZB, Delile J, Kong JH, Rousso DL, Pearson CA, et al. Olig2 and Hes regulatory dynamics during motor neuron differentiation revealed by single cell transcriptomics. PLoS Biol. 2018;16(2):e2003127.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  11. Maciel R, Bis DM, Rebelo AP, Saghira C, Zuchner S, Saporta MA. The human motor neuron axonal transcriptome is enriched for transcripts related to mitochondrial function and microtubule-based axonal transport. Exp Neurol. 2018;307:155–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Burke EE, Chenoweth JG, Shin JH, Collado-Torres L, Kim SK, Micali N, et al. Dissecting transcriptomic signatures of neuronal differentiation and maturation using iPSCs. Nat Commun. 2020;11(1):462.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Courtot AM, Magniez A, Oudrhiri N, Feraud O, Bacci J, Gobbo E, et al. Morphological analysis of human induced pluripotent stem cells during induced differentiation and reverse programming. Biores Open Access. 2014;3(5):206–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Bott CJ, Johnson CG, Yap CC, Dwyer ND, Litwa KA, Winckler B. Nestin in immature embryonic neurons affects axon growth cone morphology and Semaphorin3a sensitivity. Mol Biol Cell. 2019;30(10):1214–29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Sunabori T, Tokunaga A, Nagai T, Sawamoto K, Okabe M, Miyawaki A, et al. Cell-cycle-specific nestin expression coordinates with morphological changes in embryonic cortical neural progenitors. J Cell Sci. 2008;121(Pt 8):1204–12.

    Article  CAS  PubMed  Google Scholar 

  16. Arber S, Han B, Mendelsohn M, Smith M, Jessell TM, Sockanathan S. Requirement for the homeobox gene Hb9 in the consolidation of motor neuron identity. Neuron. 1999;23(4):659–74.

    Article  CAS  PubMed  Google Scholar 

  17. Harada A, Teng J, Takei Y, Oguchi K, Hirokawa N. MAP2 is required for dendrite elongation, PKA anchoring in dendrites, and proper PKA signal transduction. J Cell Biol. 2002;158(3):541–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Soltani MH, Pichardo R, Song Z, Sangha N, Camacho F, Satyamoorthy K, Sangueza OP, Setaluri V. Microtubule-associated protein 2, a marker of neuronal differentiation, induces mitotic defects, inhibits growth of melanoma cells, and predicts metastatic potential of cutaneous melanoma. Am J Patho. 2005;166(6):1841–50.

    Article  CAS  Google Scholar 

  19. Huang DW, Sherman BTL, R.A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4(1):44–57.

    Article  CAS  Google Scholar 

  20. Li PECP. OPaver: Omics Pathway Viewer.

  21. Tecalco-Cruz AC, Sosa-Garrocho M, Vazquez-Victorio G, Ortiz-Garcia L, Dominguez-Huttinger E, Macias-Silva M. Transforming growth factor-beta/SMAD target gene SKIL is negatively regulated by the transcriptional cofactor complex SNON-SMAD4. J Biol Chem. 2012;287(32):26764–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Huang P, Xiong F, Megason SG, Schier AF. Attenuation of notch and hedgehog signaling is required for fate specification in the spinal cord. PLoS Genet. 2012;8(6):e1002762.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Fukunaga K, Miyamoto E. Role of MAP kinase in neurons. Mol Neurobiol. 1998;16(1):79–95.

    Article  CAS  PubMed  Google Scholar 

  24. Lacinova L. Voltage-dependent calcium channels. Gen Physiol Biophys. 2005;24(Suppl 1):1–78.

    CAS  PubMed  Google Scholar 

  25. Brini M, Cali T, Ottolini D, Carafoli E. Neuronal calcium signaling: function and dysfunction. Cell Mol Life Sci. 2014;71(15):2787–814.

    Article  CAS  PubMed  Google Scholar 

  26. Tan W, Dean M, Law AJ. Molecular cloning and characterization of the human ErbB4 gene: identification of novel splice isoforms in the developing and adult brain. PLoS One. 2010;5(9):e12924.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  27. Chung DW, Wills ZP, Fish KN, Lewis DA. Developmental pruning of excitatory synaptic inputs to parvalbumin interneurons in monkey prefrontal cortex. Proc Natl Acad Sci U S A. 2017;114(4):E629–E37.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Murakami M, Ohba T, Xu F, Satoh E, Miyoshi I, Suzuki T, et al. Modified sympathetic nerve system activity with overexpression of the voltage-dependent calcium channel beta3 subunit. J Biol Chem. 2008;283(36):24554–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Takahashi K, Tanaka A, Hara M, Nakanishi S. The primary structure and gene organization of human substance P and neuromedin K receptors. Eur J Biochem. 1992;204(3):1025–33.

    Article  CAS  PubMed  Google Scholar 

  30. Steinhoff MS, von Mentzer B, Geppetti P, Pothoulakis C, Bunnett NW. Tachykinins and their receptors: contributions to physiological control and the mechanisms of disease. Physiol Rev. 2014;94(1):265–301.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Watson LM, Bamber E, Schnekenberg RP, Williams J, Bettencourt C, Lickiss J, et al. Dominant mutations in GRM1 cause Spinocerebellar Ataxia type 44. Am J Hum Genet. 2017;101(5):866.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Scott GD, Fryer AD. Role of parasympathetic nerves and muscarinic receptors in allergy and asthma. Chem Immunol Allergy. 2012;98:48–69.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Arikawa K, Takuwa N, Yamaguchi H, Sugimoto N, Kitayama J, Nagawa H, et al. Ligand-dependent inhibition of B16 melanoma cell migration and invasion via endogenous S1P2 G protein-coupled receptor. Requirement of inhibition of cellular RAC activity. J Biol Chem. 2003;278(35):32841–51.

    Article  CAS  PubMed  Google Scholar 

  34. Lukomska A, Dobrzanski G, Liguz-Lecznar M, Kossut M. Somatostatin receptors (SSTR1-5) on inhibitory interneurons in the barrel cortex. Brain Struct Funct. 2020;225(1):387–401.

    Article  CAS  PubMed  Google Scholar 

  35. Santiago C, Labrador JP, Bashaw GJ. The homeodomain transcription factor Hb9 controls axon guidance in Drosophila through the regulation of Robo receptors. Cell Rep. 2014;7(1):153–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Moore AR, Filipovic R, Mo Z, Rasband MN, Zecevic N, Antic SD. Electrical excitability of early neurons in the human cerebral cortex during the second trimester of gestation. Cereb Cortex. 2009;19(8):1795–805.

    Article  PubMed  Google Scholar 

  37. Moore AR, Zhou WL, Jakovcevski I, Zecevic N, Antic SD. Spontaneous electrical activity in the human fetal cortex in vitro. J Neurosci. 2011;31(7):2391–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Birnbaumer L, Zhu X, Jiang M, Boulay G, Peyton M, Vannier B, et al. On the molecular basis and regulation of cellular capacitative calcium entry: roles for Trp proteins. Proc Natl Acad Sci. 1996;93(26):15195–202.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Belkacemi A, Hui X, Wardas B, Laschke MW, Wissenbach U, Menger MD, et al. IP3 receptor-dependent cytoplasmic Ca(2+) signals are tightly controlled by Cavbeta3. Cell Rep. 2018;22(5):1339–49.

    Article  CAS  PubMed  Google Scholar 

  40. Decaillot FM, Kazmi MA, Lin Y, Ray-Saha S, Sakmar TP, Sachdev P. CXCR7/CXCR4 heterodimer constitutively recruits beta-arrestin to enhance cell migration. J Biol Chem. 2011;286(37):32188–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Chapman NA, Dupre DJ, Rainey JK. The apelin receptor: physiology, pathology, cell signalling, and ligand modulation of a peptide-activated class A GPCR. Biochem Cell Biol. 2014;92(6):431–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Brighton PJ, Szekeres PG, Willars GB. Neuromedin U and its receptors: structure, function, and physiological roles. Pharmacol Rev. 2004;56(2):231–48.

    Article  CAS  PubMed  Google Scholar 

  43. Lepski G, Jannes C, Nikkhah G, Bischofberger J. cAMP promotes the differentiation of neural progenitor cells in vitro via modulation of voltage-gated calcium channels. Front Cell Neurosci. 2013;7(155):1.

  44. Bando Y, Irie K, Shimomura T, Umeshima H, Kushida Y, Kengaku M, et al. Control of spontaneous Ca2+ transients is critical for neuronal maturation in the developing Neocortex. Cereb Cortex. 2016;26(1):106–17.

    Article  PubMed  Google Scholar 

  45. Dessaud E, McMahon AP, Briscoe J. Pattern formation in the vertebrate neural tube: a sonic hedgehog morphogen-regulated transcriptional network. Development. 2008;135(15):2489–503.

    Article  CAS  PubMed  Google Scholar 

  46. Nam H, Jeon S, An H, Yoo J, Lee HJ, Lee SK, et al. Critical roles of ARHGAP36 as a signal transduction mediator of Shh pathway in lateral motor columnar specification. Elife. 2019;8:e46683.

  47. Nakamura F, Kumeta K, Hida T, Isono T, Nakayama Y, Kuramata-Matsuoka E, et al. Amino- and carboxyl-terminal domains of Filamin-A interact with CRMP1 to mediate Sema3A signalling. Nat Commun. 2014;5:5325.

    Article  CAS  PubMed  Google Scholar 

  48. Dimos JT, Rodolfa KT, Niakan KK, Weisenthal LM, Mitsumoto H, Chung W, Croft GF, Saphier G, Leibel R, Goland R, Wichterle H, Henderson CE, Eggan K. Induced pluripotent stem cells generated from patients with ALS can be differentiated into motor neurons. Science. 2008;321(5893):1218–21.

    Article  CAS  PubMed  Google Scholar 

  49. McCorry LK. Physiology of the autonomic nervous system. Am J Pharm Educ. 2007;71(4):78.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Blankenship AG, Feller MB. Mechanisms underlying spontaneous patterned activity in developing neural circuits. Nat Rev Neurosci. 2010;11(1):18–29.

    Article  CAS  PubMed  Google Scholar 

  51. McGleenon BM, Dynan KB, Passmore AP. Acetylcholinesterase inhibitors in Alzheimer's disease. Br J Clin Pharmacol. 1999;48(4):471–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Pagano G, Rengo G, Pasqualetti G, Femminella GD, Monzani F, Ferrara N, et al. Cholinesterase inhibitors for Parkinson's disease: a systematic review and meta-analysis. J Neurol Neurosurg Psychiatry. 2015;86(7):767–73.

    Article  PubMed  Google Scholar 

  53. Hampel H, Mesulam MM, Cuello AC, Farlow MR, Giacobini E, Grossberg GT, et al. The cholinergic system in the pathophysiology and treatment of Alzheimer's disease. Brain. 2018;141(7):1917–33.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Kaur S, Singh S, Chahal KS, Prakash A. Potential pharmacological strategies for the improved treatment of organophosphate-induced neurotoxicity. Can J Physiol Pharmacol. 2014;92(11):893–911.

    Article  CAS  PubMed  Google Scholar 

  55. Colovic MB, Krstic DZ, Lazarevic-Pasti TD, Bondzic AM, Vasic VM. Acetylcholinesterase inhibitors: pharmacology and toxicology. Curr Neuropharmacol. 2013;11(3):315–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Lo CC, Chain PS. Rapid evaluation and quality control of next generation sequencing data with FaQCs. BMC Bioinformatics. 2014;15:366.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Shakya MC, P. PiReT Pipeline for Reference based Transcriptomics.

  58. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  59. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


We greatly valued helpful discussions with Paul Li and Bin Hu, as well as Arasely Rodriguez’s technical assistance.


Research presented in this article was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20190167.

Author information

Authors and Affiliations



ES: Methodology, Investigation, Visualization, Formal analysis, Writing – review and editing; KDA: Writing original draft, Visualization, Formal analysis; BH: Formal analysis; SMV: Formal analysis, Writing-original draft; JFH: Supervision; ST: Writing – review and editing, Supervision, and RI: Supervision, Writing – review and editing. The final version was approved by all authors.

Corresponding authors

Correspondence to Sofiya Micheva-Viteva or Rashi Iyer.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

List of differentially expressed genes (DEG) at each stage of MN development that were found to be significantly (p < 0.001) up or downregulated compared to the transcriptome profile of iPSC (D0). Additionally, DEG between two consecutive cell types are included.

Additional file 2.

Snapshot of the Ontology Pathway at the Glutamatergic synapse showing regulatory genes in cAMP and Ca2+ signaling pathways with DGE during the transition from iPSC to NSC, MNP (D13 vs D07), early MNs (D18 vs D13) and mature MNs (D28 vs D13). From left to right the color-coded legend shows DGE (p < 0.001) upon transition from iPSC to MNs. OPaver software was used to link the DGE profiles to the metabolic pathways.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Solomon, E., Davis-Anderson, K., Hovde, B. et al. Global transcriptome profile of the developmental principles of in vitro iPSC-to-motor neuron differentiation. BMC Mol and Cell Biol 22, 13 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: