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

Cell cultures from embryonic stems cells to spinal motor neurons

H9 (WA09, NIH registry 0046) hESC lines were obtained from WiCell Research Institute (Madison, WI). Human ESCs were maintained on irradiated mouse embryonic fibroblasts and differentiated as described before (Du et al., 2015). Briefly, ESCs were cultured using MN differentiation medium containing DMEM/F12, Neurobasal medium, N2, B27, ascorbic acid, Glutamax and penicillin/streptomycin (All from Gibco). The medium was additionally supplemented with chemical compounds: 3 µM CHIR99021 (Torcris), 2 µM DMH-1 (Torcris) and 2 µM SB431542 (Stemgent) for 6 days differentiated into neuroepithelial progenitors; 1 µM CHIR99021, 2 µM DMH-1, and 2 µM SB431542, 0.1 µM RA (Sigma) and 0.5 µM Purmorphamine (Pur, Sigma) for 6 days differentiated into MN progenitors; 0.5 µM RA and 0.1 µM Purmorphamine for 6 days differentiated into MNs.

Total RNA isolation and RNA sequencing

Total RNA from samples was extracted using the Trizol reagent (Invitrogen). RNA quantity and purity of total were confirmed using NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific, USA). RNA integrity was assessed by Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Ribosomal RNA was removed. Sequencing libraries were constructed. The dTTP were replaced by dUTP in the reaction buffer during second strand cDNA synthesis. Products were purified and library quality was assessed on the Agilent Bioanalyzer 2100 system. RNA Sequencing was performed on the Illumina HiSeq ((Illumina, USA) by KangChen Biotechnology Corporation (Shanghai, China) using next-generation sequencing analysis. Raw data of the RNA sequencing have been deposited in the Gene Expression Omnibus (GEO) public database (accession number: GSE151744).

Immunocytochemistry

The primary antibodies information was used: SOX1 (gIgG, 1:1000, R&D), HOXA3 (mIgG, 1:1000, R&D), OLIG2 (rIgG, 1:1000, Millipore), HB9 (rIgG, 1:1000, Millipore), NF-200 (mIgG, 1:400, CST). Image J was used to perform cell counting (NIH, USA).

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. The target siRNA sequences were as followed: 5′-ACAACCCGATGACAGCAGA-3′, 5′-CCAAATGGAGAACGTGCAA-3′and 5′-GTACTCGGTCTTTGCTGGC-3′. The target antisense oligonucleotides sequence were as followed: 5′-ATGAAAGGAAAGGCACCAGC-3′, 5′-ACATCTAACAAGGAGGACAC-3′, 5′-AGGTTGACCGCAATGCACAT-3′. The negative control (NC) Smart Silencer did not contain domains sequences homologous to those of humans, rats, or mice. The cells were transfected with the lncRNA Smart silencer using Lipofectamine 3000 (Invitrogen, USA) and collected 48 h after transfection for RNA isolation.

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 (FC) > 2.

A total of 441 DElncRNAs (192 up-regulated and 249 down-regulated), 443 DElncRNAs (198 up-regulated and 245 down-regulated) and 338 DElncRNAs (164 up-regulated and 174 down-regulated) were identified in D6 vs D0 (neuroepithelial progenitors, NEP), D12 vs D0 (motor neuron progenitors, MNP), and D18 vs D0 (motor neuron, MN) respectively (File S1). The volcano plots of expression profile for all the detected transcripts showed the relationship between the fold change and the significance (Figs. 2A–2C). To determine the key RNAs in human MN differentiation, we analyzed DElncRNAs in Venn diagram form. Veen analysis revealed that 33 lncRNAs were simultaneously up regulated and 78 were down regulated on D6, D12 and D18 (Figs. 2D–2E). Meanwhile, a total of 1068 DEmRNAs (360 up-regulated and 708 down-regulated), 1175 mRNAs (444 up-regulated and 731 down-regulated) and 68 DEmRNAs (34 up-regulated and 34 down-regulated) were identified in D6 vs D0, D12 vs D0 and D18 vs D0 respectively (File S2). The volcano plots was shown in Figs. 3A–3C. Veen analysis revealed that 8 mRNAs were simultaneously up-regulated and 34 down-regulated (Figs. 3D–3E).

Validation of lncRNAs and mRNAs expression

To verify the results of the RNA sequencing, five strongly expressed DElncRNAs and DEmRNAs during MN differentiation with FC > 40 were selected for qRT-PCR validation. As shown in Figs. 4A–4E, the expression of ZSCAN10, OCT4 and VRTN were down regulated, while HOXA3 and SP9 were up regulated. Moreover, the expression of ENST00000454596 (lncRNA rp11.001), ENST00000419695 (lncRNA rp11.003) and ENST00000583521 (lncRNA rp11-c9-001) were down regulated, while NCAM1-AS were up regulated (Figs. 4F–4J). However, lncRNA H19 increased at D12, but decreased significantly at D18. Our results agreed with the data of RNA-sequencing generally.

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 (MF) and cellular component (CC) shown in Fig. 5A. In the BP domain, the top 6 GO terms were enriched in cell migration involved in gastrulation, somatic stem 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). Bioinformatics analysis of same putative target miRNAs of HOXA3 and lncRNA NCAM1-AS identified a direct binding site of has-miR-338-3p. Thus, lncRNA-miRNA-mRNA pathway was constructed including: lncRNA NCAM1-AS-miR-338-3p-HOXA3. Moreover, 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.