Temporal expression profiles of lncRNA and mRNA in human embryonic stem cell-derived motor neurons during differentiation

Background Human embryonic stem cells (hESC) have been an invaluable research tool to study motor neuron development and disorders. However, transcriptional regulation of multiple temporal stages from ESCs to spinal motor neurons (MNs) has not yet been fully elucidated. Thus, the goals of this study were to profile the time-course expression patterns of lncRNAs during MN differentiation of ESCs and to clarify the potential mechanisms of the lncRNAs that are related to MN differentiation. Methods We utilized our previous protocol which can harvest motor neuron in more than 90% purity from hESCs. Then, differentially expressed lncRNAs (DElncRNAs) and mRNAs (DEmRNAs) during MN differentiation were identified through RNA sequencing. Bioinformatic analyses were performed to assess potential biological functions of genes. We also performed qRT-PCR to validate the DElncRNAs and DEmRNAs. Results A total of 441 lncRNAs and 1,068 mRNAs at day 6, 443 and 1,175 at day 12, and 338 lncRNAs and 68 mRNAs at day 18 were differentially expressed compared with day 0. Bioinformatic analyses identified that several key regulatory genes including POU5F1, TDGF1, SOX17, LEFTY2 and ZSCAN10, which involved in the regulation of embryonic development. We also predicted 283 target genes of DElncRNAs, in which 6 mRNAs were differentially expressed. Significant fold changes in lncRNAs (NCAM1-AS) and mRNAs (HOXA3) were confirmed by qRT-PCR. Then, through predicted overlapped miRNA verification, we constructed a lncRNA NCAM1-AS-miRNA-HOXA3 network.


INTRODUCTION
Generation of specific cell types from human embryonic stem cells (hESC) in vitro have provided powerful platforms to study human disease and to understand fundamental biological processes. Highly efficient directed differentiation of hESCs into spinal motor neurons have been used to explore not only spinal motor neuron (MN) development but also MN disorder mechanisms (Chen et al., 2014;Du et al., 2015). Motor neurons are responsible for innervating skeletal muscles in the periphery and controlling movement. Transcription factors (TFs) regulate precise temporal and spatial gene expression in motor neuron specification and differentiation (Cave & Sockanathan, 2018;Alaynick, Jessell & Pfaff, 2011). The TFs Olig2 and Ngn2 function in opposition to regulate gene expression in MN progenitors in the pMN domain and the TFs Isl1 and Lhx3 are crucial for specifying MN identity (Thaler et al., 2002;Seo, Lee & Lee, 2015;Lee et al., 2005).
LncRNAs, ranging in length from 200 nt to 100 kb, are highly expressed in the central nervous system. Accumulating evidence suggested that lncRNA played crucial roles in numerous biological and pathological processes at the chromatin remodeling level, transcriptional level and post-transcriptional level (Mercer, Dinger & Mattick, 2009;Kopp & Mendell, 2018). Notably, lncRNAs function as key regulators of cell differentiation and development, especially in neurogenesis. In particular, lncRNAs can regulate ESC pluripotency and control multiple lineage differentiation by association with miRNAs, RNA-binding proteins, and epigenetic modifiers (Loewer et al., 2010). LncRNA-1604 functioned as competing endogenous RNAs (ceRNAs) of miR-200c and indirectly regulated the core TFs ZEB1 and ZEB2 during neural differentiation from mouse ESCs (Weng et al., 2018). LncRNA Haunt functions as a genetic enhancer and an epigenetic repressor of HOXA gene activation during ESC differentiation (Yin et al., 2015). Dlk1-Dio3 locusderived lncRNAs play a critical role in maintaining postmitotic MN cell fate by repressing progenitor genes and they shape MN subtype identity by regulating Hox genes (Yen et al., 2018). Nevertheless, at present very little functional characterization of lncRNAs in human motor neuron differentiation has been elucidated.
We used spinal motor neuron differentiation to profile the temporal changes without further purification steps. We combined highly efficient MN differentiation of hESCs in vitro with RNA-seq analysis to reveal the expression profiles of lncRNAs and mRNAs. Our findings may provide a new theoretical basis for further studies on lncRNAs modulation of motor neuron differentiation.

Quantitative real-time PCR
Total RNA was extracted using the Trizol reagent (Invitrogen). qRT-PCR was then performed and the 2 − Ct method was calculated for quantification. The GAPDH was used as an internal control. The primer sequences used are listed in Table 1.

Transfection of lncRNA Smart Silencer
Motor neuron progenitors were transfected with lncRNA Smart Silencer (RiboBio Co., Guangzhou, China) to knock down the expression of NCM1-AS. The lncRNA Smart Silencer contained a mixture of three siRNAs and three antisense oligonucleotides.

Bioinformatic analysis
Gene Ontology (GO) analysis was used to investigate differentially expressed mRNAs with GO categories. The predicted target genes above were conducted using the DAVID database (http://david.abcc.ncifcrf.gov/). GO terms with a P value <0.05 were considered as significantly enriched. PPI networks was used STRING database (https://string-db.org/) and Cytoscape. The networks were visualized in CytoHubb plug-in of Cytoscape. LncRNAs and mRNAs possessing microRNA recognition elements (MREs) for the targeted miRNAs were predicted using the miRanda and TargetScan.

Statistical analysis
All qRT-PCR results are expressed as the means ± SEM of at least three independent experiments. Statistical analyses were performed with SPSS statistics software version 22.0. P values < 0.05 was considered statistically significant. All graphs were made with GraphPad Prism 8.

Differentiation of high purity motor neuron from human embryonic stem cells
In this study, we used hESCs to differentiate into spinal cord MNs in vitro. The MNs were generated using chemical protocol described as method part (Du et al., 2015) (Fig. 1A). The hESCs can differentiate into SOX1 neuroepithelial progenitors at day 6, OLIG2 motor neuron progenitors at day 12, and HB9 motor neuron at day 18 high efficiently by using a combination of small molecules (Figs. 1B-1D).

Patterns of gene expression changes from hESCs to motor neuron
High throughput sequencing is an efficient approach for investigating the biological function of RNAs. All the differently expressed DElncRNAs and DElncRNAs were statistically significant (P < 0.05) with fold change (

Bioinformatics analysis during MN differentiation
To identify the key factors that regulated MN differentiation, GO analysis of DEmRNAs was performed on three different aspects namely biological process (BP), molecular function ( cell population maintenance, positive regulation of transcription from RNA polymerase II promoter, positive regulation of cell proliferation, cardiac cell fate determination and anterior/posterior pattern specification. Genes associated with cell migration involved in gastrulation were SOX17, MIXL1 and CER1. Genes associated with positive regulation of cell proliferation were EPHA1, ETS1, HOXA3, POU3F3, FLT1, and TDGF1. In the CC domain, the top 3 GO terms were transcription factor complex, membrane raft and nucleoplasm. In the MF domain, the top 3 GO terms were sequence-specific DNA binding, transcription factor activity and HMG box domain binding. Figure 5B showed the heat map of DEmRNAs. Furthermore, the Protein-Protein Interaction (PPI) network of DEmRNAs contained 28 nodes and 39 edges (Fig. 5C). The topological analysis of the network was carried out by Network Analyzer in Cytoscape. PPI network was imported into Cytohubba to determine the hub transcription factors with high degree of connectivity between the nodes. The top ten hub genes were POU5F1, TDGF1, SOX17, LEFTY2, ZSCAN10, CER1, ZFP42, MIXL1, L1TD1 and ESRP1 shown in Fig. 5C.

Analysis of lncRNAs target mRNAs
LncRNA may regulate nearby protein-coding genes by cis-regulatory effects. We analyzed potential function of the target genes of DElncRNAs at D18 vs D0. Total 338 DElncRNAs had 289 mRNAs target genes (File S3). Then we performed GO analysis on DElncRNAs target genes to explore their potential biological functions. The GO terms related to BP, CC and MF were shown in Fig. 6A. The most enriched BP term was related to positive regulation of telomerase activity. As for CC, we found that the most enriched term were nucleoplasm and transcription factor complex. As for MF, transcription factor activity, sequence-specific DNA binding and metal ion binding were most enriched. Moreover, KEGG analysis was made in Fig. 6B and the most enriched 5 top pathway terms were MAPK signaling pathway, Ribosome, Pyrimidine metabolism, Hippo signaling pathway and Cushing syndrome. Veen diagram analysis indicated that 6 mRNAs targeted by DElncRNAs were found at the interaction of DEmRNAs at D18 vs D0 (Fig. 6C). Moreover, lncRNAs and their cis target DEmRNAs were up-regulated at D18, including two well-known TFs HOXA6 and HOXC9. We also performed the topological analysis of PPI network on these target genes (Fig. 7), indicating the highly connected hub nodes in PPI network.

Construction of lncRNA-miRNA-mRNA interaction network
LncRNAs can regulate gene expression by acting as ceRNAs to sponge miRNAs (Thomson & Dinger, 2016). Thererfore, we constructed a ceRNA interaction network from verified DE mRNAs and DE lncRNAs based on previous qRT-PCR data. As the RNA sequencing and PCR data shown, the transcription factors HOXA3 and SP9, exhibited continuous up-regulation in the transition from ESCs to MN stages, and especially showed a sharp up-regulation coincident with MNs specification. Interestingly, the expression of lncRNA NCAM1-AS showed same trend as HOXA3 and SP9. The NCAM1-AS is antisense lncRNA with two exons, produced from gene NCAM1, which is involved in cell adhesion, axonal outgrowth, synapse formation during development and differentiation, and highly expressed in the developing central and peripheral nervous systems (Wobst et al., 2015). the expression of HOXA3 was down-regulated at D12 in the NCAM1-AS Smart Silencer group compared to NC group (Fig. 8A). Furthermore, immunofluorescence analysis showed a decreased expression of OLIG2 in the NCAM1-AS knockdown group (Fig. 8B). Therefore, the results suggested that lncRNA NCAM1-AS can affect MNP differentiation.

DISCUSSION
Motor neuron differentiation is precisely regulated and orchestrated by combinatorial expression of TFs during embryogenesis (Alaynick, Jessell & Pfaff, 2011). Accumulating evidence suggested that lncRNAs could interact with transcription factors to regulate cell differentiation (Lopez-Pajares et al., 2015;Ng et al., 2013;Wang et al., 2014). In this study, we profiled mRNAs and lncRNAs expression from our highly efficient ESC-derived MN differentiation protocol to study the development of MNs. Our analysis focused on the identification of transcription factors and lncRNAs that are strongly involved in the temporal development of MNs. ESCs, derived from the inner cell mass of blastocyst stage embryos, can both self-renew and differentiate into other cell types (Thomson et al., 1998). The balance between selfrenewal and differentiation is regulated by a complex interaction network of translation factors. Pluripotent genes such as Oct4, Nanog, Sox2, Klf4, and Myc (Chambers & Tomlinson, 2009), activated in ESCs, were inhibited during cell differentiation, whereas expression of differentiation marker genes increases gradually.
Notably, we found the well-known pluripotency-associated transcription factors POU5F1 (also known as OCT4) and TDGF1 which were hub downstream-regulated genes upon MN differentiation in our study. It has been reported that lincRNA linc-RoR, may function as a key ceRNA to link the network of miRNAs and core TFs OCT4, SPX2, and NANOG, thus regulating ESC maintenance and differentiation (Wang et al., 2013). As a previous study (Zhang et al., 2006), hub transcription factor ZSCAN10, verified by qRT-PCR, also down-regulated in our study and could regulate ESCs gene expression and differentiation. OCT4 can directly regulate expression of ZSCAN10 and TDGF1 (Wang et al., 2007;Babaie et al., 2007).
Here, we also identified cis regulatory target genes HOXA6 and HOXC9 of lncRNA HOXA-AS3 and HOXC-AS2 at MN stage, respectively. Spinal MNs acquire specialized pool identities that guide their axons to target muscles in the limb, and the specificity of these precise connections (Dasen et al., 2005). MNs could express many HOX genes specifying MN pool identity and connectivity (Dasen et al., 2005;Lacombe et al., 2013). HOX6 paralog group genes (HOXA6, HOXC6, and HOXB6) contributed to diverse aspects of motor neuron subtype differentiation, and determined lateral motor column (LMC) fate at forelimb levels of the spinal cord (Lacombe et al., 2013). In addition, HOXC9 determined thoracic level MN population fates, including preganglionic column (PGC) and hypaxial motor column (HMC) neurons (Jung et al., 2010). The LncRNA HOXA-AS3 was found to inhibit osteogenic differentiation and promote adipogenic differentiation (Wu et al., 2017). The up-regulated lncRNAs HOXA-AS3 and HOXC-AS2 at human MN derived from ESCs was first identified in our study. The two lncRNAs might be involved in MN differentiation by cis-regulating HOXA6 and HOXC9, which needs further study.
Our sequencing results suggested a series of mRNAs and lncRNAs significantly changed during the transition from ESC to motor neurons. HOXA3, SP9 and lncRNA NCAM1-AS verified by PCR were observed dramatically up-regulated, especially at period of motor neuron. HOXA3 and HOXB3 are necessary for the specification of Pax6-and Olig2dependent somatic MN progenitors (Gaufo, Thomas & Capecchi, 2003). In addition, HOX1 was reported to be involved in mediating both the role of RA-signaling in specification of hindbrain MNs (Schubert et al., 2006). In our sequencing data, HOXB3 was up-regulated at NEP and MNP stages but not altered at MN stage. The neuronal differentiation marker NCAM was involved in motor neurons functionally expanding synaptic territory (Chipman, Schachner & Rafuse, 2014).
Additionally, the potential function of lncRNA NCAM1-AS has been originally identified in human MN differentiation from ESC. The lncRNAs could function as miRNA sponges and might compete against other endogenous RNAs to regulate mRNA expression levels and maintain normal biological function. Bioinformatics analysis indicated lncRNA-miRNA-mRNA pathway: lncRNA NCAM1-AS-miR-338-3p-HOXA3. A recent study found that miR-338-3p targeted and inhibited HOXA3 in breast cancer (Zhang & Ding, 2019). Thus, the upregulated lncRNA NCAM1-AS might inhibit the expression of miRNA by acting as a miRNA sponge, and in turn, increasing the expressions of MN differentiation-associated mRNAs HOXA3. The dysregulated lncRNAs may regulate gene expression through many ways and play a critical role in the processes of neuronal differentiation.

CONCLUSIONS
In conclusion, we utilized our highly efficient ESC-derived MNs differentiation protocol and next-generation sequencing to provide new insights into understanding the molecular mechanisms underlying MN differentiation and modulating lineage commitment of ESCs. The understanding of MN differentiation could ultimately offer the early diagnosis and