Single-cell RNA-seq reveals dynamic transcriptome profiling in human early neural differentiation

Abstract Background Investigating cell fate decision and subpopulation specification in the context of the neural lineage is fundamental to understanding neurogenesis and neurodegenerative diseases. The differentiation process of neural-tube-like rosettes in vitro is representative of neural tube structures, which are composed of radially organized, columnar epithelial cells and give rise to functional neural cells. However, the underlying regulatory network of cell fate commitment during early neural differentiation remains elusive. Results In this study, we investigated the genome-wide transcriptome profile of single cells from six consecutive reprogramming and neural differentiation time points and identified cellular subpopulations present at each differentiation stage. Based on the inferred reconstructed trajectory and the characteristics of subpopulations contributing the most toward commitment to the central nervous system lineage at each stage during differentiation, we identified putative novel transcription factors in regulating neural differentiation. In addition, we dissected the dynamics of chromatin accessibility at the neural differentiation stages and revealed active cis-regulatory elements for transcription factors known to have a key role in neural differentiation as well as for those that we suggest are also involved. Further, communication network analysis demonstrated that cellular interactions most frequently occurred in the embryoid body stage and that each cell subpopulation possessed a distinctive spectrum of ligands and receptors associated with neural differentiation that could reflect the identity of each subpopulation. Conclusions Our study provides a comprehensive and integrative study of the transcriptomics and epigenetics of human early neural differentiation, which paves the way for a deeper understanding of the regulatory mechanisms driving the differentiation of the neural lineage.

Background: Investigating cell fate decision and subpopulation specification in the context of the neural lineage is fundamental to understanding neurogenesis and neurodegenerative diseases. The differentiation process of neural-tube-like rosettes in vitro is representative of neural tube structures, which are composed of radially organized, columnar epithelial cells and give rise to functional neural cells. However, the underlying regulatory network of cell fate commitment during early neural differentiation remains elusive. Results: In this study, we investigated the genome-wide transcriptome profile of single cells from six consecutive reprogramming and neural differentiation time points and identified cellular subpopulations present at each differentiation stage. Based on the inferred reconstructed trajectory and the characteristics of subpopulations contributing the most towards commitment to the central nervous system (CNS) lineage at each stage during differentiation, we identified putative novel transcription factors in regulating neural differentiation. In addition, we dissected the dynamics of chromatin accessibility at the neural differentiation stages and revealed active cis-regulatory elements for transcription factors known to have a key role in neural differentiation as well as for those that we suggest are also involved. Further, communication network analysis demonstrated that cellular interactions most frequently occurred among embryoid body (EB) stage and each cell subpopulation possessed a distinctive spectrum of ligands and receptors associated with neural differentiation which could reflect the identity of each subpopulation. Conclusions: Our study provides a comprehensive and integrative study of the transcriptomics and epigenetics of human early neural differentiation, which paves the way for a deeper understanding of the regulatory mechanisms driving the differentiation of the neural lineage.

89
The advance of single cell trans-omics technology has offered incisive tools for 90 revealing heterogeneous cellular contexts and developmental processes [9][10][11].

154
It is widely reported that chromatin structures undergo widespread reprogramming 155 during cell status transition, with some genomic regions becoming compacted or 156 opened, leading to the switching on or off of a repertoire of genes responsible for cell 157 fate decision [24][25][26][27][28][29]. We studied the dynamic chromatin landscape by tracing the 158 temporal origins of ATAC peaks at each stage with peaks non-overlapping with 159 existing ones that were annotated as novel peaks. We assumed that those peaks, 160 conserved among differentiation stages, are associated with housekeeping genes 161 while stage-dynamic peaks are likely to represent cis-regulatory elements important 162 for cell status transition. As expected, we observed the introduction of roughly 10-50% 163 of novel peaks in each stage, accompanied by the disappearance of several 164 pre-existing ATAC peaks. Notably, more novel peaks appeared at the NPCs stage 165 than at other stage ( Fig. 1c). GO term analysis of genes residing in novel peaks 166 across the differentiation stages showed enrichment of "axon development", "positive 167 regulation of nervous system development", "epithelial tube morphogenesis", "positive 168 regulation of neurogenesis", "cell-cell signalling by Wnt", "forebrain development", 169 "hindbrain development", "telencephalon development", "neural precursor cell 170 proliferation", and "cell fate commitment". "Neurotrophin signalling pathway" was also 171 found to be enriched, but was specifically associated with NPCs. KEGG enrichment 172 analysis showed that "FoxO signalling pathway", a pathway which is known to play an 173 important role in NPC proliferation, and "neuroactive ligand−receptor interaction" were enriched in NPCs stage (Fig. 1d,e), suggesting that specific cis-regulatory elements 175 regulating neural differentiation are being staged (poised) for stem cell fate 176 specification and conversion.

178
To reveal the detail of chromatin accessibility dynamics during neural differentiation, 179 we also analysed the gained or lost peaks at each stage compared with the previously 180 neighbouring one. We observed that the number of gained peaks was with the largest 181 increase at the NPCs stage while the number of lost peaks was relatively high at 182 Ros-E stage (Additional file 2: Figure S2a). Next, we studied the genomic distribution 183 of these dynamic peaks and found that both the gained and lost peaks were located 184 mostly in distal intergenic regions and promoter regions (Additional file 2: Figure S2b).

185
This observation indicates that distal and promoter regions are more dynamic 186 compared to other genomic regions during neural differentiation process.

188
To gain insight into the potential function of closing (lost) peaks dynamics, we carried 189 out GO enrichment analysis on the genes associated with lost peaks at each stage.

329
The ectoderm marker, OTX1, and genes involved in the ventral hindbrain marker (e.g.,

469
Among the 131 TFs exhibiting differential expression from Ros-L3 to NPC1, 80 TFs 470 were up-regulated while 51 TFs were down-regulated (Additional file 11: Figure S11; 471 Additional file 19:  Figure S11). This is consistent with our previous observation that the main 478 trajectory has progressed more towards to CNS. Of particular interest is PRDM1, 479 whose expression increased from Ros-E2 to Ros-L3 and decreased during the 480 progression from Ros-L3 to NPC1 ( Fig. 4g and Additional file 10, 11: Figure S10, 11),

481
suggesting that it might play multiple specific roles in neural differentiation. To dissect the cis-regulatory elements directing the expression of those regulators, we 495 selected the differentially expressed TFs that showed differential ATAC peaks 496 between neighbouring stages and performed motif scanning on the differential peaks.

497
Focusing on the transition from Ros-E2 to Ros-L3, we found transcription factor 498 binding sites (TFBSs) for TEAD2 and YY1 in a differential ATAC peak downstream of 499 the PRDM1 gene (Fig. 4c). Multiple motifs for the transcription factor TFAP2C were 500 found in a differential peak located in the intron of the ARID3A gene, which is a 501 regulator responsible for the transition for Ros-L3 to NPCs (Fig. 4d). Based on the 502 temporal specificity of ATAC peaks and the existence of TF motifs in these regions, we 503 propose that those elements are stage-specific cis-regulatory elements regulating the 504 expression of neural regulators in response to their upstream regulatory TFs.

506
To infer the putative targets of key regulators, we combined the information from ATAC 507 peaks and motifs for TFs. All peaks containing motifs for a certain TF were annotated 508 as TF-related peaks and genes proximal to the peak were considered as potential 509 targets of that TF. Using these criteria, we predicted thousands of targets for each TF 510 (Additional file 20: Table S2). To dissect the regulatory network of each TF, we 511 conducted GO term and KEGG enrichment analysis for the putative target list of each 512 key regulator. Our results suggested that, from Ros-E2 to Ros-L3, the targets for 513 PRDM1 were significantly enriched in pathways and GO terms associated with "axon 514 guidance", "hippo signalling pathway" and "neurotrophin signalling pathway" (Fig. 4e   515 and Additional file 13: Figure S13a). From Ros-L3 to NPC1, targets for NR2F1, SOX9, 516 TFAP2C were enriched in KEGG pathways associated with "axon guidance" and 517 "hippo signalling pathway" (Additional file 13: Figure S13b- these cell lines (Fig. 4g, h).  Table S3). The most frequent interactions were observed in the EB stage, implying 534 that cells communicate extensively to coordinate differentiation programs during 535 embryogenesis (Additional file 14: Figure S14). In contrast, much fewer interactions   can be used to confidently infer cell identity.

573
During development in vivo the neuroectoderm folds to form the neural tube which is 574 then patterned into regionally specialized subunits composed of progenitor cells.

749
Together, our findings are supported by different genetic cell lines mitigating the 750 concern that our results are limited to the cells forming the basis of this study.

752
Through differential expression analysis, we identified genes specifically expressed at 753 each stage which include both cell status master regulators such as TFs and  Neural differentiation 804 We applied a well-adopted neural differentiation protocol [8,16] Assay for transposase-accessible chromatin sequencing (ATAC-seq) 846 We profiled open chromatin accessibility sequencing (ATAC-seq) of neural

862
The filtered data were aligned to the reference genome (hg19)    Additional file 1: Figure S1. Quality control of ATAC-seq.
Color scheme is based on z-score distribution from -2 (purple) to 2 (yellow). Gene symbols highlight with color specific to the respective Ros-L subset. b Box plots of selected TFs defined in Figure S7a. c Top 12 of GO terms identified by up-regulated genes specific to the respective Ros-L subpopulation with the color as indicated (red: GO terms specific to Ros-L1; purple: GO terms specific to Ros-L2; blue: GO terms specific to Ros-L3; gray: selected GO terms shared by Ros-L1 and Ros-L3). d KEGG enrichment analysis of Ros-L2 (all terms) and Ros-L3 (selected terms), respectively.
-1 0 1 2 Z-score Additional file 8: Figure S8 a NPC2 NPC3 of discriminative TF sets for each cluster in NPCs stage with P-value cutoff ≤ 0.01.
Color scheme is based on z-score distribution from -1 (purple) to 2 (yellow). Gene symbols highlight with color specific to the respective NPC subset. b Box plot of selected TFs defined in Figure S8a. c Top 10 (NPC1) and all (NPC2 and NPC3) of GO terms identified by up-regulated genes specific to the respective Ros-L subpopulation with the color as indicated (blue: GO terms specific to NPC1; green: GO terms specific to NPC2; pink: GO terms specific to NPC3).   iPSCs to EB3 EB3 to Ros-E2 Additional file 12: Figure S12. Key regulators during neural differentiation. a

Regulatory network of differentially expressed TFs between iPSCs and EB3. b
Regulatory network of differentially expressed TFs between EB3 and Ros-E2.