A single-cell transcriptional roadmap for cardiopharyngeal fate diversification

Abstract

In vertebrates, multipotent progenitors located in the pharyngeal mesoderm form cardiomyocytes and branchiomeric head muscles, but the dynamic gene expression programmes and mechanisms underlying cardiopharyngeal multipotency and heart versus head muscle fate choices remain elusive. Here, we used single-cell genomics in the simple chordate model Ciona to reconstruct developmental trajectories forming first and second heart lineages and pharyngeal muscle precursors and characterize the molecular underpinnings of cardiopharyngeal fate choices. We show that FGF-MAPK signalling maintains multipotency and promotes the pharyngeal muscle fate, whereas signal termination permits the deployment of a pan-cardiac programme, shared by the first and second heart lineages, to define heart identity. In the second heart lineage, a Tbx1/10-Dach pathway actively suppresses the first heart lineage programme, conditioning later cell diversity in the beating heart. Finally, cross-species comparisons between Ciona and the mouse evoke the deep evolutionary origins of cardiopharyngeal networks in chordates.

Figure 1∣. Cell clustering and cell-type-specific markers.

(a) Early cardiopharyngeal development in Ciona, and sampling stages and established lineage tree. Cardiopharyngeal lineage cells are shown for only one side and known cell-type-specific marker genes are indicated. st., FABA stage; hpf, hours post-fertilization; TVC, trunk ventral cell; STVC, second trunk ventral cell; FHP, first heart precursor; ASMF, atrial siphon muscle founder cells; SHP, second heart precursor; iASMP, inner atrial siphon muscle precursor; oASMP, outer atrial siphon muscle precursor; LoM, longitudinal muscles; QC, quality control. Dotted line: midline. (b) t-distributed Stochastic Neighbor Embedding (t-SNE) plots of 20 hpf scRNA-seq data (n=288 cells) showing distinct clusters of progenitor subtypes: FHP (red), SHP (orange), iASMP (blue) and oASMP (dark blue). Color-coded marker gene expression levels are shown on corresponding clusters. (c) Expression heatmap of 20 hpf single cell transcriptomes showing top predicted differentially expressed marker genes across different cell types. Blue: previously known ASM and heart markers, red: candidate markers. (d) Violin plots and FISH validations of candidate cell-type-specific markers in St. 28 embryos. mRNAs visualized by whole mount fluorescent in situ hybridization (green). Cardiopharyngeal nuclei marked by Mesp>nls::LacZ revealed by anti beta-galactosidase antibody (red). Mesp>hCD4::mCherry, revealed by anti-mCherry antibody, marks cell membranes (blue). Anterior to the left. Scale bar, 10 μm. Solid arrowheads, ASM; open arrowheads, SHPs; arrows, FHPs; M, midline (dotted line). The numbers of observed embryos and those showing the illustrated gene expression pattern are indicated at the right bottom corner of each image. Violin plots are to visualize the distributions of the expression (log FPKM) of the indicated genes. The wide of the violin indicates the frequency of cells with indicated gene expression level. The number of cells in each cell cluster is summarized in Supplementary Table 6 (Data sheet: cell identity and number).

Figure 2∣. Reconstruction of cardiopharyngeal developmental trajectories.

(a) Cell lineages used to reconstruct three unidirectional cardiopharyngeal trajectories. (b) Diffusion maps showing the cardiopharyngeal trajectories. Color-coded cell identities as defined by unsupervised clustering from larvae dissociated at indicated time points (Supplementary Fig. 1a). Black lines: principal curve; light gray contours: single cell density distribution. Color codes correspond to assigned cell identities following clustering at each time point. hpf, hours post-fertilization. DC: Diffusion Coordinate. (c) Distribution of identified cell types isolated at defined time points along the trajectories, showing the general agreement between the time series and developmental progression, but also that cells isolated from a given time point are not all at the same developmental “pseudotime”. (d) Cross-correlation heatmaps to infer regulatory states along the trajectories. Dendrogram (left) obtained from constrained hierarchical clustering. Top bars indicate the sample of origin with color codes as in (c). PCC, Pearson Correlation Coefficient. (e) Relative cell identity composition for each regulatory states identified on the trajectories. Note the 16ASM cells clustering with the ‘STVC’ state in the ASM trajectory, indicating that these cells retain most STVC characteristics and have not yet activated the ASM-specific program. (f) Pseudotemporal expression profiles of indicated genes along the FHP trajectory. X-axis: normalized pseudotime as defined in (b), Y-axes: relative expression level. Black lines indicate the smoothed expression. Black dashed lines indicate the transitions between predicted regulatory states as defined in (d) and color-coded below. Purple dashed lines indicate calculated activation or inactivation pseudotime. Dot colors refer to the sample of origin as indicated in (c). (g) Proportions of primed vs. de novo-expressed genes among defined categories of marker genes.

Figure 3∣. Transcriptional regulation of ASM fate specification.

(a) Pseudotemporal expression profiles of indicated genes along the ASM trajectory. X-axis: normalized pseudotime as defined in Fig.2. Y-axes: relative expression levels. Black lines indicate the smoothed expression. Red dashed lines indicate the transitions between predicted regulatory states, indicated as in Fig. 2d, and purple dashed lines indicate calculated activation pseudotime. (b) Heatmap of smoothed expression profiles along the ASM trajectory showing primed pan-cardiac genes (inhibited, left) and candidate de novo ASM genes (activated, right). White vertical lines mark transitions between indicated regulatory states along the ASM trajectory. Colored bars on the left indicate the regulatory state of calculated activation pseudotime. (c) Violin plots showing the log₂(fold-change) of candidate ASM-specific genes (n=159) in response to indicated perturbations of Ebf function, a dominant-negative (Ebf^DN) and Ebf over-expression (Ebf^OE) as in Razy-Krajka et al.. The white bars indicate the interquartile range. The black whiskers extended from the bars represent the upper (max) and lower (min) adjacent values in the data. The black lines in the middle of the bars show the median values. (d) Chord diagram showing mutual enrichment of ASM vs. Cardiac genes among candidate target genes activated or inhibited by Ebf, respectively. Ebf is predicted to downregulate a few ASM candidate genes, which are primed and quickly downregulated after ASM specification (e.g. Hand-r⁴). (e) Heatmap of smoothed expression profiles for candidate ASM-specific Ebf target genes defined in Razy-Krajka et al., showing activation pseudotimes in regulatory states ASM2 and ASM3 (left blue bars). White vertical lines mark transitions between indicated regulatory states. (f) Predicted induction pseudotime of candidate ASM-specific Ebf target genes (black dashed lines separate corresponding ASM regulatory states).

Figure 4∣. A pan-cardiac program for heart fate specification.

(a) Smoothed gene expression along FHP pseudotime for de novo-expressed pan-cardiac genes (activated) and primed ASM genes (down-regulated). White vertical lines: transitions between predicted regulatory states. (b) Principal component (PC1) correlates with pseudotime. PCC: Pearson’s Correlation Coefficient. Sample size (cells on the FHP trajectory) n=379. (c) Average PC1-loading scores for indicated gene category, mapped onto the FHP trajectory. (d) Proportions of de novo pan-cardiac and FHP-specific genes with calculated activation pseudotime in binned pseudotime windows along FHP trajectory. (e) Expression profiles of Mmp21 along the FHP trajectory. Purple dashed line: calculated activation pseudotime. Dots colors: samples of origin as in Fig. 2b. (f) Activation pseudotimes for de novo-expressed pan-cardiac genes along first and second heart lineages pseudotime axes. (g) Proportions of de novo-expressed pan-cardiac genes with calculated activation pseudotime in binned pseudotime windows along FHP and SHP trajectories. (d-e, g) X-axis: normalized pseudotime. Black dashed lines: transitions between regulatory states.

Figure 5∣. FGF-MAPK signaling regulates the pan-cardiac program for heart fate specification.

(a) Differential expression of primed and de novo-expressed ASM and pan-cardiac genes in indicated conditions vs. LacZ. (b) Violin plots represent the distributions of log₂ fold changes in indicated conditions and time points relative to LacZ controls, and parsed by primed or de novo-expressed pan-cardiac genes. The white bars indicate the interquartile range. The black whiskers extended from the bars represent the upper (max) and lower (min) adjacent values in the data. The black lines in the middle of the bars show the median values. Sample size: Primed Pan-cardiac genes, n= 58; de novo Pan-cardiac genes, n= 26. Summary statistics: results of two-tailed t-test for significant difference from 0 are indicated below violin plots, results for KS tests for significant differences between “Primed” and “de novo” gene sets in each condition are indicated above violin plots. “=”, no difference, “**”, P-value <0.01. (n = 2 biological replicates; Supplementary Table 6). (c) FGF-MAPK inhibition induces precocious Lrp4/8 expression in multipotent cardiopharyngeal progenitors (TVCs). Lrp4/8 mRNAs (green) visualized by FISH and processed by Imaris (green dots). Anti-beta-galactosidase antibody (red) marks TVC nuclei expressing Mesp>nls::LacZ. Mesp>hCD4::mCherry, revealed by anti-mCherry antibody (blue), marks cell membranes. Anterior to the left. Scale bar, 10 μm. Box plots represent the distributions of numbers of Lrp4/8+ dots per cell in indicated conditions. Bars indicate the median value. * p=4.46e⁻¹¹ ( One-tailed student’s t-test, n = 2 biological replicates). (d) Summary model showing the maintenance and progressive restriction of FGF-MAPK signaling in the multipotent progenitors (TVCs, trunk ventral cells, and STVCs, second trunk ventral cells) and atrial siphon muscle founder cells (ASMFs). Inhibition of MAPK activity permit the deployment of de novo-expressed pan-cardiac genes in both cardiac lineages (FHP, first heart precursors, and SHP, second heart precursors). FHPs specifically activate genes like Mmp21, and later produce most Mhc2+ cardiomyocytes, whereas SHPs descend from Tbx1/10+ multipotent progenitors, and activate Dach, which contributes to inhibiting the FHP-specific program. See discussion and Razy-Krajka et al. for details.)

Figure 6∣. A second-heart-lineage specific Tbx1/10-Dach pathway.

(a) Smoothed gene expression along the SHP trajectory for SHP-specific genes. White vertical lines: transitions between regulatory states. (b) Dach expression pattern along the SHP trajectory. Purple dashed line: predicted induction time. Black dashed lines: transitions between regulatory states. (c) Perturbations of Dach function does not alter Ebf expression in ASMPs. Numbers: observed/total. (d) Dach and Tbx1/10 are required to restrict Mmp21 activation to the FHP. Barplots: proportions of larvae showing indicated phenotypes in each experimental condition. Wt, wild-type. Foxf>Dach^PAMmut: TVC-specific Foxf enhancer driving expression of CRISPR-resistant Dach cDNA. (e) Tbx1/10 function is required for Dach expression in SHP. Panels 2 and 4: segmented Dach+ dots superimposed on cell patterns. (c-e) Confocal stacks acquired for 10 larvae in each condition in biological duplicates. None of the 20 Tbx1/10^CRISPR larvae showed Dach expression in SHPs. Solid arrowheads: ASMPs, open arrowheads: SHPs, arrows: FHPs. (f) FGF-MAPK signaling negatively regulates Dach expression in Tbx1/10+ ASMFs. Representative confocal stacks showing segmented Dach+ dots in 18.5 hpf (St. 27) larvae. Blocking MEK activity causes ectopic Dach expression in the ASMFs (solid arrowheads), in addition to its endogenous expression in the SHPs (open arrowheads) (n=10/10 for each condition). (g) Combined Tbx1/10 over-expression and MAPK inhibition induced precocious Dach expression in 14 hpf B7.5 lineage cells. Open arrows: STVCs, arrows: FHPs; dotted line: midline (m). Violin plots represent the distribution of counts Dach+ dots per cell. Black dots: cells with identify and experimental perturbation indicated below. Red dots: mean values. The thin red line: the upper (max) and lower (min) adjacent values in the data. One-tailed student’s t-test indicates the precocious Dach expression in both FHPs and STVCs in combined Tbx1/10 over-expression and MAPK inhibition condition. * p=9.972e-05; # p=0.005275. (c-g) mRNAs (green) visualized by FISH or segmented (green dots). Mesp>nls::LacZ, revealed by anti-beta-galactosidase antibody (red), marks nuclei. Mesp>hCD4::mCherry, revealed by anti-mCherry antibody (blue), marks cell membranes. Dotted line: midline (m). Anterior to the left. Scale bar, 10 μm.

Figure 7∣. Characteristics and origins of the intracardiac cell diversity in beating hearts.

(a) Lineage tracing by photoconversion of nuclear Kaede. Mesp>nls::Kaede::nls (green) marks nuclei of live B7.5 lineage cells. Kaede photoconverted from green (I) to red (II) specifically in the FHPs of a 16hpf larva, which is shown at successive time points (III and IV). Segmented nuclei are shown below. Open arrowheads, ASMPs; solid arrowheads, SHPs; arrows, FHPs, Dotted line: midline. arrowheads, SHPs; arrows FHPs; Anterior to the left. Scale bar, 10 μm. Experiment performed in biological duplicates. (b) Confocal data (left) and segmented image (right) showing Mhc2 expression (grey) primarily excluded from SHP-derived cells in wild-type juvenile heart (St. 38). Mesp>nls::LacZ marks nuclei of FHP and SHP-derived cells (red), Tbx1/10>H2B::mCherry marks only SHP-derived cells (green). Scale bar, 10μm. (c) Dach and Tbx1/10 antagonize the production of Mhc2+ cardiomyocytes from the second heart lineage. Rendered segmented signals are shown. Grey: Mhc2 mRNA. Green: Tbx1/10>H2B::mCherry, revealed by an anti-mCherry antibody, marks SHP-derived cells; Red: Mesp>nls::LacZ, revealed by an anti-beta-galactosidase antibody, marks all B7.5 lineage cells. Scale bar, 10μm. Boxplots: proportions of Mhc2+ cells among the Tbx1/10>H2B::mCherry+ SHP-derived cells in juvenile hearts. Bars in the box indicate the median value. * p=1.41e⁻⁴ and * p=2.80e⁻⁵ for Dach^CRISPR and Tbx1/10^CRISPR, respectively (One-tailed student’s t-test). N numbers represent the embryos analyiszed for the experimental perturbation as indicated. (d) t-SNE plots of scRNA-seq data acquired in n=386 FACS-purified cardiopharyngeal lineage cells from juveniles at St. 38. (e) Feature plots: expression of indicated markers in clusters shown in (d). (f) Top predicted differentially expressed genes across the st. 38 juvenile heart. (g) Matn1/3/4 and Dach are enriched in the Tbx1/10>H2B::mCherry+ SHP-derived cells (green). Scale bar, 5μm. Images are XY cross section of juvenile hearts. Boxplots: proportions of Matn1/3/4+ or Dach+ cells among indicated cell populations. Both Matn1/3/4 (p = 4.459e-11) and Dach (p = 0.006795) are significantly enriched in SHP-derived cells. N numbers represent juveniles used to quatifiy the gene expression among indicated cell populations. Bars in the box: median value in each condition. * p<0.05 (One-tailed student’s t-test).

Figure 8∣. Conserved cardiopharyngeal programs in chordates.

(a) t-SNE plots of mouse scRNA-seq data (n=1,205 cells) showing Tbx1 and Dach1 expression patterns. Cluster identities are as determined in the original publication, with the pharyngeal mesoderm shown in green. (b) Expression patterns of Dach1, Islet1 and Nkx2.5 proteins in E9.5 mouse embryos. Representative image of 4 analysed embryos. Arrows: Islet1+ head muscle progenitor cells in the mesodermal core of the first (PA1, top) and second (PA2, bottom) pharyngeal arches, showing absence of Dach1 and Nkx2.5 expression. Open arrowhead: Dach1+, Islet1+, Nkx2.5− second heart field cells in the dorsal pericardial wall, solid arrowhead: Triple Nkx2.5+, Dach1+, Islet1+ second heart field-derived cells in the outflow tract (OFT). Note the Nkx2.5+, Dach1−, Islet1− cells in the ventricle (V). Scale bar, 100 μm. (c) Aligned structure of Ciona and Mouse E8.25 cardiopharyngeal cells (n=2291 cells). tSNE plots showing the clustering of Ciona and Mouse E8.25 cardiopharyngeal cells, respectively, using conserved markers determined by canonical correlation (CC). Barplots indicate that original cell identities, defined in each species independently, as recovered in the clustering using conserved markers. (d) tSNE plots of Ciona and Mouse scRNA-seq data as described in (c), with the expression patterns of Ebf1 and Gata4. (e) Single cell expression profiles for the top 30 conserved markers in each species, separately.