Inferring Regulatory Programs Governing Region Specificity of Neuroepithelial Stem Cells during Early Hindbrain and Spinal Cord Development

Abstract

Neuroepithelial stem cells (NSC) from different anatomical regions of the embryonic neural tube's rostrocaudal axis can differentiate into diverse central nervous system tissues, but the transcriptional regulatory networks governing these processes are incompletely understood. Here, we measure region-specific NSC gene expression along the rostrocaudal axis in a human pluripotent stem cell model of early central nervous system development over a 72-h time course, spanning the hindbrain to cervical spinal cord. We introduce Escarole, a probabilistic clustering algorithm for non-stationary time series, and combine it with prior-based regulatory network inference to identify genes that are regulated dynamically and predict their upstream regulators. We identify known regulators of patterning and neural development, including the HOX genes, and predict a direct regulatory connection between the transcription factor POU3F2 and target gene STMN2. We demonstrate that POU3F2 is required for expression of STMN2, suggesting that this regulatory connection is important for region specificity of NSCs.

Keywords: context-specific regulatory networks; early human neural development; gene regulation; transcriptional regulatory networks.

Figure 1. Analysis of time course data for inference of context-specific dynamic regulatory networks for region-specific neuroepithelial stem cells (NSCs).

A. Shown are the major steps of our workflow. Starting with mRNA RNA-seq samples from hindbrain and cervical spinal cord NSCs (HB/CSC NSCs), Escarole was used to cluster genes into discrete expression states per time point (i), enabling discovery of global transcriptional dynamics (ii) and dynamic genes (iii). In parallel, we inferred a context-specific regulatory network (iv) from the HB/CSC NSC RNA-seq, publicly available RNA-seq data from similar neural stem cell types, and a DNase I-based network structure. We integrated the regulatory network with the dynamic genes identified by Escarole to prioritize key regulators (v). B. Overview of context-specific network inference. First, we apply MERLIN-P to learn a general neural lineage network (All-NSC) by integrating RNA-seq from this study, published studies and a DNase I-based footprint network. Second, the HB/CSC data are used again within MERLIN-P, taking the All-NSC network as a prior, to infer a HB/CSC-specific regulatory network (Refined-NSC). C. Schematic timeline of Hox propagation and NSC differentiation for RNA-seq dataset generation (STAR Methods). Gray bars depict onset and duration of changes to culture media. (i) hPSCs are differentiated to region-specific NSCs. HOX propagation begins at day 2, at which point samples are taken at intervals of 2 hours (*) from 0h-24h and 6 hours (†) from 30h-72h. (ii) Each sample is exposed to retinoic acid (RA) for four days to induce differentiation to the neural stem cell state. D. Representative image of PAX6+ neuroepithelial cells forming rosettes in vitro. Red = PAX6, Green = N-cadherin (CDH2); Scale bar = 100um. E. Principal Components Analysis of the RNA-seq time course data. Each point is labeled with its time point. F. The top five genes most positively correlated to PC 1, ranked by their coefficients (loadings).

Figure 2. Escarole modeling of time series data.

A. Probabilistic graphical model of Escarole. Left shows the graphical model for three cell types. S_A, S_B and S_C denote the state variables for the cell types, A, B and C, and y_A, y_B, and y_C, denote the observed expression variables. The transition probabilities are denoted by P(S_x|S_y) for each pair of cell types, x, y ∈ {A, B, C} and the emission probabilities are denoted by P(y_x|S_x). Right shows further details of the model with two expression states (interpreted as UP and DN) and three linearly ordered cell types. In each cell type c, a Gaussian distribution $N\left({\mu }_x^s,{\sigma }_x^s\right)$ models the emission probability of a gene’s expression value. Dynamics of state transitions, P(S_x|S_y) between cell types are captured by the white black transition matrices, which give the prior probability of a gene’s expression state given its state in the parent cell type. B. Simulation-based comparison of Escarole to related methods for identifying dynamic genes. Results are averaged over five simulated data sets. Positives are defined as genes that ever transition their expression state. ‘# genes’ gives the average total number of predicted dynamic genes. C. Ability of methods to predict timing of gene transitions (when the simulated expression state is different for the x-axis time point and the previous time point). Points represent means over five trials, error bars ± SD. D. Ranking of neural relevant gene sets when algorithms are applied to the HB/CSC RNA-seq data from this study. Height of bars is mean adjusted p-value over five runs of GSEAPreranked; ± SD. HOX: HOX transcription factors; Others: genes important for differentiation to neuroepithelial cells (NE), early radial glial (ERG), mid-radial glial (MRG), or any (UNION) (Ziller et al., 2015). E. Mean Euclidian silhouette index (SI) of genes within per-timepoint clusters. Points represent means over five trials; error bars ± SD. F. Distributions of per-gene SI for global clusters or trajectory genesets that span all time points, excluding singletons. See also Figure S1, Table S1.

Figure 3. Escarole analysis of NSC RNA-seq time course.

A. Escarole expression states for each time point/cell type of the RNA-seq time course. Each row of heat maps corresponds to the five expression states at one time point; each column shows that expression state (1–5) across all time points. The height of each heat map is proportional to the number of genes assigned to the expression state at the row time point. B. Boxplots showing the expression values from genes assigned to each state, aggregated over all time points. See also Figure S2. C, D. Gene ontology, KEGG, REACTOME and cis-regulatory element enrichment analysis of genes assigned to Escarole states within each time point. See also Table S2. Shown are the top 5 pathways (C) or cis-regulatory elements (D) with FDR-corrected hypergeometric p<1E-5, per expression state in each time point, ranked by corrected p-value. Each block of rows shows one of the five Escarole states; each column shows one time point.

Figure 4. Analysis of Escarole state trajectories.

See also Table S3. A. Expression values and Escarole inferred states for HOXB9. Expression in heatmaps is mean-centered per gene for visualization. The line plot shows the uncentered values. B. Bird’s eye view of 53 trajectory gene sets composed of genes with similar trajectories defined by hierarchical clustering. The right-hand columns indicate whether a trajectory gene set is significantly enriched for members of a curated gene set (e.g., Gene Ontology, KEGG Pathways), targets of a transcription factor, or a gene set targeted by chemical toxins. Bolded module IDs are referenced in the text. C. Gene expression profiles and Escarole states of genes in set #1405. Columns with purple entries to the left indicate which genes are targets of a regulator; columns to the right indicate membership in curated pathways and interactions with neurotoxins for which the gene set is enriched. Bottom heatmaps show enriched regulators with available mRNA levels in the NSC RNA-seq time course. D. Members of gene set #1517; legend is the same as for C.

Figure 5. Dynamic genes identified by Escarole analysis.

A. Shown are 92 genes that transition at least 2 expression states and are highly expressed (state 4 or 5) at some point. Gene names in red are mentioned in the text. See also Table S4. B. Connections between 67 dynamic genes and eight neurotoxins whose gene targets are enriched for dynamic genes (FDR-corrected p<0.05).

Figure 6. Inferred context-specific transcriptional regulatory network for hindbrain, spinal cord NSCs.

See also Table S5, Dataset S1. A. Visualization of the Refined-NSC network. Node sizes are proportional to out-degree (number of outgoing connections), with the top 11 regulators labeled. B. Out-degree distribution of regulators. C. Assessment of inferred regulatory targets for individual TFs by comparison to transcription factor ChIP-seq or knockdown in mouse or human neural cell types. AUPR: Area under the Precision-Recall Curve. EP: targets inferred based on enhancer promoter interactions; DE: differential expression; MRG: Mid-radial glial cell type; NE: neuroepithelial cell type. D. GSEAPreranked enrichment of regulators for important neural differentiation genes identified by shRNA screens. Regulators were ranked by their importance for predicting expression of 4,931 genes that transition Escarole expression state. (i) Performance of regulator ranking by the Refined-NSC network. (ii) Performance of regulator ranking by the All-NSC-NPC network. E. Inferred interactions amongst HOX genes in the Refined-NSC network. Color/arrowhead shape indicates the sign of the expression correlation between the two genes (adjusted p<0.05).

Figure 7. Regulator prioritization and experimental validation.

A. Regulators from Refined-NSC network prioritized for dynamic genes, with inset showing the top 22 regulators (rank of maximum enrichment for HOX genes (red text)). See also Table S6. B. Subnetwork of Refined-NSC network connecting dynamic genes to their inferred regulators. Edge color and arrow heads denote activating or repressive relations. C. Subnetworks from B showing the two-neighborhoods of HOXA5 (i, #1), MAF (ii, #2), and POU3F2 (iii, #5). Right: Heatmaps show the gene expression and Escarole states of the regulators and their dynamic target genes. Additional columns mark transcription factors (‘TF’) and show the scores of the edge from the Refined-NSC, All-NSC network, and NPC motif network (‘NPC PIQ’). D. Schematic of RNAi knockdown experiments (STAR Methods). Cells are transitioned to retinoic acid at Day 4 to differentiate to neural stem cells, matching the 48hr time point in the RNA-seq experiment (Figure 1C). E. mRNA levels (qRT-PCR) for candidate regulators, targets, and control genes following treatment with non-targeting siRNA (siNT, dark blue) or RNAi-mediated knockdown of GAPDH (siGAPDH, aqua), HOXA (siHOXA5, green), or POU3F2 (siPOU3F2, yellow). x-axis labels indicate the measured gene. Error bars ± SD. Statistical analysis was performed using a one-way ANOVA with Tukey’s posthoc analysis (n=3 biological replicates per condition). *p< 0.05; **p< 0.01. See also Table S7.