Single-cell analysis of bidirectional reprogramming between early embryonic states identify mechanisms of differential lineage plasticities in mice

Abstract

Two distinct lineages, pluripotent epiblast (EPI) and primitive (extra-embryonic) endoderm (PrE), arise from common inner cell mass (ICM) progenitors in mammalian embryos. To study how these sister identities are forged, we leveraged mouse embryonic stem (ES) cells and extra-embryonic endoderm (XEN) stem cells-in vitro counterparts of the EPI and PrE. Bidirectional reprogramming between ES and XEN coupled with single-cell RNA and ATAC-seq analyses showed distinct rates, efficiencies, and trajectories of state conversions, identifying drivers and roadblocks of reciprocal conversions. While GATA4-mediated ES-to-iXEN conversion was rapid and nearly deterministic, OCT4-, KLF4-, and SOX2-induced XEN-to-induced pluripotent stem (iPS) reprogramming progressed with diminished efficiency and kinetics. A dominant PrE transcriptional program, safeguarded by GATA4, alongside elevated chromatin accessibility and reduced DNA methylation of the EPI underscored the differential plasticities of the two states. Mapping in vitro to embryo trajectories tracked reprogramming cells in either direction along EPI and PrE in vivo states, without transitioning through the ICM.

Keywords: ES cells; XEN cells; blastocyst; epiblast; extra-embryonic endoderm; lineage plasticity; pluripotency; primitive endoderm; reprogramming; single-cell analysis.

Figure 1.. Oct4, Sox2 and Klf4 can successfully reprogram XEN cells in a slow and inefficient manner.

(A) XEN cells were derived from mice homozygous for 3 alleles (R26:M2rtTA^Tg/Tg; Col1a1:OKSmCh^Tg/Tg; Oct4:EGFP^Tg/Tg referred to as ‘3F’). (B) (Left) Experimental scheme describing the reprogramming conditions used for XEN cells. (Right) Brightfield and fluorescent images of XEN cells during the reprogramming time course. Scale bars represent 250μm. (C) RT-qPCR data for several XEN and ES markers in wild-type XEN (IM8A-1), 3F XEN, two XEN-iPS lines, and wild-type ES cells (R1). Individual bars show mean expression of three technical replicates normalized to mean expression of two reference genes: Actb and Gapdh; error bars represent standard deviation. (D) Unsupervised hierarchical clustering of bulk RNA-seq data from XEN-iPS, wild-type ES, and XEN cells. (E) (Top) Experimental scheme outlining the generation of XEN-iPS chimeric embryos. (Bottom) Maximum intensity projection of 5 optical sections from a confocal image of a chimeric embryo stained with indicated markers. Nuclei labeled with DAPI, mCherry counter stained with anti-RFP. Scale bars represent 100μm.

Figure 2.. Reciprocal lineage conversions of XEN and ES cells have drastically different kinetics and efficiencies of conversion.

(A) Time course flow cytometry-based analysis of pluripotency-associated markers SSEA-1 and Oct4-GFP, and XEN-associated marker PDGFRα during XEN reprogramming. Representative contour plots show expression of SSEA-1 (AlexaFluor647-conjugated), Oct4-GFP and PDGFRα (PE-Cy7-conjugated) at days 2, 8 and 16 of reprogramming. Population percentage indicated within each gate. (B) (Top) Tracking the reprogramming efficiencies of four major subpopulations sorted at day 14 of reprogramming. At day 28, populations arising from each sorted subpopulation were determined using flow cytometry. (Bottom, left) Stacked bar charts depicting the mean proportion of the population represented by each subpopulation at day 14, and at day 28 (bottom, right). S-O-P+ subpopulation constitutes the remainder of the population (not shown) to amount to 100% and represents non-reprogramming XEN cells. N = 6 (2 independent experiments with 3 replicates each). (C) Time course tracking of ES-to-iXEN conversion using flow cytometry analysis of SSEA-1, PDGFRα and Gata6-Venus reporter. Representative contour plots show expression of PDGFRα (PE-Cy7-conjugated), Gata6-Venus and SSEA-1 (BV421-conjugated) at indicated timepoints. Population percentage is indicated within each gate. (D) Summarized schematic of XEN-to-iPS and ES-to-iXEN lineage conversions. (Left) Hypothesized reprogramming route taken by XEN cells. Sizes and number of arrows reflect the likelihood of cells progressing from one state to the next (smaller/single arrows = few cells progress; larger/more arrows = many cells progress). Reprogramming XEN cells initiate expression of Oct4-GFP, then downregulate PDGFRα. This downregulation step represents a bottleneck during reprogramming since few Oct4-GFP+ cells progress to this state. Following PDGFRα downregulation, a large proportion of cells upregulate SSEA-1. (Right) Hypothesized lineage conversion of ES to iXEN cells. Initial upregulation of both, PDGFRα and Gata6-Venus, is followed by downregulation of SSEA-1. No obvious bottlenecks are detected during ES-to-iXEN conversion.

Figure 3.. scRNA-seq analyses of XEN-to-iPS and ES-to-iXEN conversions.

(A) Experimental scheme outlining cells and lineage conversion timepoints assayed by scRNA-seq. (B) Force-directed layouts showing XEN-to-iPS (top) and ES-to-iXEN (bottom) conversion trajectories. Individual plots highlight individual timepoints; pie charts indicate the proportion of the trajectory represented by each timepoint. (C) Palantir determined pseudotime ordering, differentiation potential, and terminal states of ES-to-iXEN (top) and XEN-to-iPS (bottom) trajectories. (D) Branch probabilities of terminal states determined by Palantir in the XEN-to-iPS trajectory. Black arrowhead indicates cells at T1 with low probability of differentiating to T2. Green arrowhead indicates where T2 probability increases. Red arrowhead indicates cells at T2, having a non-zero probability of acquiring the T3 state. (E) Gene expression of XEN and pluripotency-associated markers. Each cell is colored by its MAGIC imputed expression level. T1, T2 and T3 terminal states indicated for the XEN-to-iPS trajectory.

Figure 4.. XEN-to-iPS and ES-to-iXEN conversion trajectories approximate in vivo cell states.

(A) Force-directed layout of combined XEN-to-iPS and ES-to-iXEN trajectories based on Harmony integration (number of cells = 91,447). XEN-to-iPS (left) or ES-to-iXEN (right) are highlighted. Cells are colored by Palantir pseudotime, computed separately as in Figure 3. (B) Gene expression of pluripotency and XEN-associated markers displayed in the combined trajectory. Each cell is colored by its MAGIC imputed expression level in the individual XEN-to-iPS or ES-to-iXEN trajectories, as labeled. (C) Force-directed layout of combined trajectories (as in panel A) with individual timepoints from each trajectory colored as indicated. (D) (Top) in vivo embryo stages and tissues (labeled) profiled by scRNA-seq. (Bottom) Force-directed layouts of combined in vivo stages and lineages as labeled above, and color-coded as indicated by cell type. (E) Gene expression signatures of XEN, XEN-iPS, and T1 and T2 terminals states mapped onto force-directed layouts of combined in vivo states including (top) or excluding (bottom) the EPI lineage. (F) Visualization of in vitro XEN-to-iPS (left) or ES-to-iXEN (right) bins mapped onto the combined in vivo trajectory. Individual circles represent single in vitro bins and are color coded according to pseudotime.

Figure 5.. XEN transcriptional network presents a roadblock for XEN-to-iPS reprogramming.

(A) Gene expression trends over pseudotime for XEN-to-iPS reprogramming. Plots show mean expression trend of all genes within each cluster. Dotted curve indicates probability of acquiring the T3/iPS state. Vertical dotted lines indicate T1 and T2 terminal states along the pseudotime axis. Representative genes for each cluster are highlighted in boxes. N = number of genes in each cluster. (B) Reprogramming efficiency following Gata4 perturbation. XEN cells were transfected with plasmids expressing Cas9 and sgRNAs targeting Gata4 or control. Immunofluorescence staining determined successful knockout of GATA4 expression (panel C). Following 10 days of reprogramming, resulting percentage of iPS-like cells was determined using flow cytometry (panel D). (C) Box plots showing relative reduction in anti-GATA4 fluorescence immunostaining in transfected XEN cells. Individual points represent relative fluorescence intensity in individual cells as arbitrary units and normalized to untransfected (RFP-) cells within the same well (see images in supplementary Figure S5C). (D) Box plots showing proportion of iPS-like cells (S+O+P-) in the entire population at day 10 of reprogramming following GATAT4 KO. N = 3. (E) Experimental timeline of GATA4 KO and sample collection for RNA-seq 10 days after the start of reprogramming. (F) Volcano plot showing significantly up/downregulated genes in cells transfected with Gata4-targeting sgRNA compared with control. (G) Stacked bar chart depicting the percentage of downregulated, upregulated, and non-differentially expressed genes that belong to the gene expression trend clusters shown in (A). (H) Heatmap representation of gene expression changes (z-score) between empty vector and Gata4 sgRNA transfected cells (two different sgRNAs were used, G4–1 and G4–2), and association with the clusters shown in (A). (I) Tornado plots of normalized ATAC-seq signal in XEN and ES cells around cell type-specific peaks. ATAC-seq signals are shown for 2.5kb up/downstream of peak centers. (J) Box plots showing relative chromatin accessibility (normalized ATAC-seq signal) in XEN and ES cells of genes enriched in XEN, XEN-iPS, T1 and T2 terminal states.

Figure 6.. Establishing an EPI-like chromatin state underlies the inefficient XEN-to-iPS conversion.

(A) Venn diagrams depicting number and percentages of highly variable accessible peak loci identified from scATAC-seq data for XEN-to-iPS conversion trajectory that overlap with differential accessible peak loci identified from bulk ATAC-seq data for XEN and ES cells (top), or with H3K27ac ChIP-seq peak loci in XEN and ES cells (bottom). (B) Force-directed layouts of combined scATAC-seq dataset for XEN-to-iPS conversion, highlighting individual timepoints (top left), metacell clusters (top right), or ChromVAR scores for XEN- and ES-specific TFs (bottom). Number of cells = 61,040; number of metacells = 543. (C) (Top) Force-directed layout of XEN-to-iPS conversion scATAC-seq dataset highlighting the groups of metacells identified based on similar accessibility profiles along pseudotime. (Bottom) Heat map view of relative chromatin accessibility changes over pseudotime of XEN-to-iPS conversion for a subset of selected peak loci, from group 1 to group 5 of the metacells. Chromatin accessibility is measured in units of number of fragments (median-normalized) across metacells followed by log transformation with a pseudo-count of 1. Each row corresponds to a peak locus and each column corresponds to a metacell. Peak loci shown are differential accessible peak loci identified in the XEN-to-iPS conversion based on differential accessibility analysis of peak loci between group 4B and group 5 of the metacells. (D) (Top) Force-directed layout showing the scATAC-seq metacell data for ES-to-iXEN conversion and highlighting the individual timepoints (number of cells = 95,396; number of metacells = 567). (Bottom) Heatmap view of the relative chromatin accessibility changes of ES-to-iXEN conversion based on the highly variable peak loci that overlap between the ES-to-iXEN and XEN-to-iPS conversions. (E) Representative IGV (Integrative Genomics Viewer) tracks showing accessibility peaks of pseudo-bulk scATAC-seq data of XEN-to-iPS (top) or ES-to-iXEN (bottom) conversion. Highlighted are relative accessibility in metacell groups (top panel; XEN-to-iPS) or individual timepoints (bottom panel; ES-to-iXEN) at example genomic loci showing late opening in XEN-to-iPS, versus gradual opening in ES-to-iXEN conversion. Signal values are indicated to the right. (F) LOLA enrichment analysis of late opening, gradual opening, late closing or transient opening peaks during XEN-to-iPS conversion. (G) Boxplots showing differences in DNA methylation levels between XEN and ES cells of highly variable loci depicted in (C). (H) XEN reprogramming with different doses of 5-Azacytidine. Box plots show the proportion of iPS-like cells (S+O+P-) at day 10 for reprogramming. N = 3.