Gene regulatory network reconfiguration in direct lineage reprogramming

Abstract

In direct lineage conversion, transcription factor (TF) overexpression reconfigures gene regulatory networks (GRNs) to reprogram cell identity. We previously developed CellOracle, a computational method to infer GRNs from single-cell transcriptome and epigenome data. Using inferred GRNs, CellOracle simulates gene expression changes in response to TF perturbation, enabling in silico interrogation of network reconfiguration. Here, we combine CellOracle analysis with lineage tracing of fibroblast to induced endoderm progenitor (iEP) conversion, a prototypical direct reprogramming paradigm. By linking early network state to reprogramming outcome, we reveal distinct network configurations underlying successful and failed fate conversion. Via in silico simulation of TF perturbation, we identify new factors to coax cells into successfully converting their identity, uncovering a central role for the AP-1 subunit Fos with the Hippo signaling effector, Yap1. Together, these results demonstrate the efficacy of CellOracle to infer and interpret cell-type-specific GRN configurations, providing new mechanistic insights into lineage reprogramming.

Keywords: cell fate prediction; direct lineage reprogramming; gene perturbation simulation; gene regulatory networks; machine learning; single-cell analysis.

Graphical abstract

Figure 1

Application of CellOracle to assess reprogramming GRN dynamics (A and B) Overview of CellOracle. (A) First, CellOracle uses scATAC-seq data to identify accessible regulatory elements, which are scanned for TF binding motifs, generating a Base GRN—a list of potential regulatory connections between a TF and its target genes (B). (C) Using single-cell expression data, active connections are identified from all potential connections in the base GRN. (D) Cell type- and state-specific GRN configurations are constructed by pruning insignificant or weak connections. (E) Hnf4α and Foxa1-mediated fibroblast to iEP reprogramming. (F) (Left) Force-directed graph: 15 clusters of cells are grouped into five cell types; fibroblasts (Fib), early transition (Early), transition (Tran), dead-end, and reprogrammed iEPs (iEP). (Right) Projection of Apoa1 (iEP marker) and Col1a2 (fibroblast marker) expression. (G) CellOracle analysis. Heatmap (left) and boxplot (right) of network edge strength between Hnf4α-Foxa1 and its target genes. ∗∗∗p < 0.001. (H) Degree and eigenvector centrality scores for Hnf4α-Foxa1. (I) Hnf4α-Foxa1 network cartography terms for each cluster. (J and K) Scatterplots of degree centrality scores between specific clusters. (J) Degree centrality score comparison between Fib_1 cluster GRN and other early and transition reprogramming cluster GRNs. (K) Degree centrality score comparison between iEP_1 and Dead-end_0 cluster GRNs.

Figure 2

Lineage tracing links early network state to reprogramming outcome (A) Overview of CellTag-based clonal tracking. Cells are transduced with the random CellTag lentiviral library so that each cell expresses three to four CellTags, resulting in a unique, heritable barcode signature. CellTags are transcribed and captured during single-cell profiling, enabling clonally related cells to be tracked throughout an experiment. (B) Experimental strategy to capture state-fate relationships. MEFs are transduced with Hnf4α-Foxa1 for 48 h, then transduced with CellTags. The end of this period is considered reprogramming day 0. Cells are expanded, and 25% of the population is profiled at day 4; this is termed the state population. The remaining cells are reseeded and profiled again on days 10 and 28 to capture reprogramming fate. (C) Captured state-fate cells. Time point information projected onto the Uniform Manifold Approximation and Projection (UMAP) embedding. A total of 24,799 cells were sequenced: 8,440 on day 4, 4,836 on day 10, and 11,523 on day 28. (D) Projection of fibroblast, iEP, and dead-end identity scores and (E) fate annotations onto the UMAP embedding. (F) A randomized test identified day 4 state clones whose day 10 and 28 fate sisters were iEP-enriched or iEP-depleted. (Top) Kernel density estimation of iEP-enriched day 4 state clones and their day 10 and 28 fates, outlining the reprogramming trajectory (n = 1,347 cells). (Bottom) iEP-depleted state-fate cells outlining the dead-end trajectory (n = 4,802 cells). (G) Projection of iEP-enriched and iEP-depleted clones onto the UMAP embedding. (H) Comparison of degree centrality scores between native fibroblasts and day 4 reprogrammed-destined cells (left) and day 4 reprogrammed- and dead-end-destined cells (right).

Figure 3

Systematic in silico simulation of TF KO to identify novel regulators of iEP reprogramming (A) Monocle-based pseudotemporal ordering of 48,515 cells from Biddy et al. (2018), two independent biological replicates. (B) Schematic for perturbation score calculations. CellOracle calculates a perturbation score by comparing the direction of the simulated cell state transition with the direction of cell differentiation. First, the pseudotime data is summarized by grid points and converted into a 2D gradient vector field. The results of the perturbation simulation are converted into the same vector field format, and the inner product of these vectors is calculated to produce a perturbation score. (C) A positive perturbation score (green) suggests that the perturbation is predicted to promote reprogramming. In contrast, the negative perturbation score (magenta) represents impaired reprogramming. (D) Ranked list of TFs based on the sum of the negative perturbation score. (E) Representative examples of TF KO simulation (top row). Expression of respective genes (bottom row). (F) Experimental validation of candidate TFs: colony-formation assay. (G) Colony quantification. n = 5 independent biological replicates for non-targeting scramble shRNA control, Fosb, Id1; n = 4 independent biological replicates for Eno1, Klf4; n = 3 independent biological replicates for Fos; unpaired t test with Welch's correction, two-tailed; ^∗p < 0.05, ^∗∗p < 0.01.

Figure 4

CellOracle analysis and experimental validation of Fos in establishing and maintaining iEP identity (A) Degree centrality, betweenness centrality, and eigenvector centrality of Fos for each cluster. (B) Network cartography terms of Fos for each cluster. (C) Fos expression projected onto the force-directed graph. (D) Violin plot of Fos expression across reprogramming stages. ∗∗∗p < 0.001. (E and F) (E) Fos gene overexpression simulation with reprogramming GRN configurations. (Left) The projection of simulated cell transitions onto the force-directed graph. The Sankey diagram summarizes the simulation of cell transitions between cell clusters. For overexpression simulation, Fos expression was set to 1.476, representing its maximum value in the imputed gene expression matrix (F) Fos gene KO simulation. (G) Colony-formation assay with addition of Fos to Hnf4α-Foxa1. (Left) E-cadherin immunohistochemistry. (Right) Boxplot of colony numbers (n = 6 technical replicates, two independent biological replicates; ^∗∗∗p < 0.001, t test, one sided). (H) qPCR for Fos and iEP marker expression (Apoa1 and Chd1) following addition of Fos to Hnf4α-Foxa1 (n = 3 independent biological replicates; ^∗∗∗p < 0.001, ^∗∗p < 0.01, t test, one sided). (I) Fos gene KO simulation in expanded, long-term cultured iEPs. (J) CRISPR-Cas9 Fos KO in expanded iEP cells. (Left) Kernel density estimation was applied with the t-SNE (t-distributed stochastic neighbor embedding) to compare cell density between control guide RNAs and guide RNAs targeting Fos. (Right) Quantification of changes in cell ratio following Fos KO.

Figure 5

Inferred Fos targets reveal a role for the Hippo signaling effector, Yap1, in reprogramming (A) Heatmap of expression of the top 50 inferred Fos targets across reprogramming. Established YAP1 targets are highlighted in red. (B) Colony-formation assay with the addition of Yap1 and Fos to Hnf4α-Foxa1. (Left) E-cadherin immunohistochemistry. (Right) Boxplot of colony numbers (n = 6 independent biological replicates; ^∗∗∗p < 0.001, t-test, one sided). (C) Brightfield and epifluorescence images of cells reprogrammed with Hnf4α-Foxa1 or Hnf4α-Foxa1-Fos-Yap1. Scale bar, 500 μm. (D) scRNA-seq of cells reprogrammed with Hnf4α-Foxa1 (n = 7,414 cells), Hnf4α-Foxa1-Fos (n = 8,771 cells), Hnf4α-Foxa1-Yap1 (n = 8,549 cells), and Hnf4α-Foxa1-Fos-Yap1 (n = 10,507 cells), profiled at day 20. Projection of fibroblast and iEP identity scores onto the UMAP embedding. (E) Kernel density estimation of cell density for each reprogramming cocktail from (D). (F) Violin plot of iEP identity scores for each reprogramming cocktail. ^∗∗∗∗p < 0.0001, Wilcoxon test. (G) Unsupervised cell type classification for each reprogramming cocktail, using normal and injured mouse liver as a reference. BEC, biliary epithelial cells. ^∗p < 0.0001, randomized test.