An Orthologous Epigenetic Gene Expression Signature Derived from Differentiating Embryonic Stem Cells Identifies Regulators of Cardiogenesis

<div><p>Here we used predictive gene expression signatures within a multi-species framework to identify the genes that underlie cardiac cell fate decisions in differentiating embryonic stem cells. We show that the overlapping orthologous mouse and human genes are the most accurate candidate cardiogenic genes as these genes identified the most conserved developmental pathways that characterize the cardiac lineage. An RNAi-based screen of the candidate genes in <i>Drosophila</i> uncovered numerous novel cardiogenic genes. shRNA knockdown combined with transcriptome profiling of the newly-identified transcription factors zinc finger protein 503 and zinc finger E-box binding homeobox 2 and the well-known cardiac regulatory factor NK2 homeobox 5 revealed that zinc finger E-box binding homeobox 2 activates terminal differentiation genes required for cardiomyocyte structure and function whereas zinc finger protein 503 and NK2 homeobox 5 are required for specification of the cardiac lineage. We further demonstrated that an essential role of NK2 homeobox 5 and zinc finger protein 503 in specification of the cardiac lineage is the repression of gene expression programs characteristic of alternative cell fates. Collectively, these results show that orthologous gene expression signatures can be used to identify conserved cardiogenic pathways.</p></div

ZNF503, ZEB2 and NKX2-5 are necessary for differentiation of CMs from human ESCs.

<p>(A) Flow cytometric analysis of cTnT expression on human ESCs after 10 days of differentiation along the cardiac lineage following shRNA knockdown of the indicated genes. Representative results of five experiments. WT cells were infected with a shRNA targeting GFP or empty vector, which yielded identical results as non-infected cells. (B) Hierarchical clustering of genes showing significant gene expression changes following shRNA knockdown of the indicated genes at the indicated time points. WT = wild-type; ZE = ZEB2 shRNA; ZN = ZNF503 shRNA; NK = NKX2-5 shRNA.</p

A multi-species epigenetic and gene expression signature identifies candidate cardiogenic genes.

<p>(A) The distribution of the indicated histone modifications to the genomic region of the mouse and human NKX2-5 gene at the ESC or CP state. Genomic coordinates are indicated for human chromosome 5 (hg19) and mouse chromosome 17 (mm9). Brackets indicate scale of peak. (B) Percentage of indicated genes known to regulate mammalian cardiogenesis. (C) Percentage of indicated cardiogenic genes conserved in <i>Drosophila</i>. (D) Enriched GO categories associated with the indicated genes. * p < 0.0001; **p < 0.0005.</p

Schematic of enhancer classification beginning with a small training set.

<p>A small set of known enhancers active in similar cells is increased by incorporating orthologous sequences. <i>De novo</i> and known motifs are mapped onto this training set and a set of control sequences. Feature vectors are used to build an enhancer model based upon the learned motif weighting. This model can be used to scan the genome for similar enhancers as the training set. These predictions can be tested using transgenic reporter assays or analysis of the expression of the associated gene. The motif weighting can likewise be used to identify novel classes of transcriptional regulators. The role of the motifs can be tested in <i>cis</i> and the identification of co-regulating TFs can be subsequently tested in <i>trans</i>.</p

The wild-type activities of FC enhancers require input from classifier-defined Myb and POUHD TF binding motifs.

<p>(A) TRANSFAC position weight matrices for Myb (V$MYB_Q6) and POUHD (V$POU1F1_Q6) enriched motifs identified by the classifier. (B) Binding site sequences in the <i>Ndg</i> enhancer for Myb and POUHD and versions in which those sites are selectively mutated. Motifs were defined by searching for matches to the vertebrate homologues in the UniPROBE database <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Robasky1" target="_blank">[99]</a>. The identification of these binding sites and the designs of the mutant versions are described in <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531.s013" target="_blank">Table S4</a>. (C) GFP (green) and β-Gal (magenta) are co-expressed when driven by the wild-type (WT) <i>Ndg</i> enhancer (<i>Ndg^WT-GFP</i> and <i>Ndg^WT-lacZ</i>, respectively). (D) GFP (green) expression driven by a version of the <i>Ndg</i> enhancer in which POUHD sites are selectively inactivated (<i>Ndg^POUHD-GFP</i>) is significantly reduced compared to β-Gal (magenta) driven by <i>Ndg^WT-lacZ</i>. (E) β-Gal driven by a version of the <i>Ndg</i> enhancer in which Myb binding sites are selectively inactivated (<i>Ndg^Myb-lacZ</i>) is de-repressed into additional somatic mesodermal cells compared to GFP driven by a WT version of the <i>Ndg</i> enhancer (<i>Ndg^WT-GFP</i>).</p

Candidate enhancers predicted by the classifier are active in FCs.

<p><i>In situ</i> hybridization of <i>dve</i> in wild-type (WT) embryos and embryos over-expressing Ras (Twi>Ras) in the mesoderm (A). Note the increased activity of <i>dve</i> in Twi>Ras embryos, indicative of a FC gene <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Estrada1" target="_blank">[22]</a>. GFP driven by the classifier-predicted enhancers associated with the upstream sequences of <i>slou</i> (arrows in B) and <i>slp1</i> (arrows in C). Slou protein (magenta) co-expresses with GFP (green) in <i>slou-GFP</i> embryos (B). Duf (magenta), which marks all FCs, co-expresses with <i>slp1</i>-<i>GFP</i> (green) (C). GFP (D) driven by the classifier-predicted intronic sequence associated with the <i>dve</i> gene co-expresses with Mef2 (D′) in myotubes at stage 15 in <i>dve-GFP</i> embryos.</p

DNA binding domains of the TFs most relevant to FC enhancer classification.

<p>Only DNA binding domains for the fifty most relevant TFs have been included. TFs were ranked according to the SVM weights of their respective motifs, which represent their discriminating power. Only the highest scoring motif for each TF was considered (median ranks computed across 10 random partitions of the training data varied between 12 and 117). <i>De novo</i> motifs were explicitly excluded from this analysis. TF domains and sequences have been clustered using average linkage and Euclidean distance. The dendogram on top of the heatmap represents the relationships among the sequences in the training data, built on the presence/absence of TFBSs recognized by a specific class of TF DNA binding domain. The dendogram on the left of the heatmap shows the relationships among the different TF DNA binding domains.</p

The enhancer classifier performs with high specificity and sensitivity.

<p>(A) Over-representation of TFBSs in the training set including only <i>D. melanogaster</i> enhancers and in the set extended using phylogenetic profiling, as compared with background sequence. P-values were adjusted for multiple testing using the method of Benjamini and Hochberg (BH) <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Benjamini1" target="_blank">[120]</a>. (B) Average ROC curve for the 10-fold cross-validation. Our method achieves an area under the ROC curve of 0.89 (shaded in gray). FPR: false-positive rate; TPR: true-positive rate. (C) Distribution of FC enhancer scores for the genome-wide scan. Scores assigned by the classifier for each evaluated sequence are shown in red. We used a FPR of 5% to define a cut-off for putative enhancers (dotted blue line; see <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#s4" target="_blank">Materials and Methods</a> for details). (D) Fold-enrichment in 180 validated FC genes in the neighborhood of putative FC enhancers, as determined for different FPRs. Intergenic putative FC enhancers were associated with the closest gene, whereas intronic sequences were associated with their host gene. P-values were computed using the binomial test.</p

TFBS combinatorics within FC enhancers.

<p>(A) Distribution of Tcf, Mad, Pnt, Twi, Tin, POUHD, Tbx, Myb, Fkh, HD and Mef2 TFBSs in FC enhancers. Binding sites for Tcf, Mad, Pnt, Twi and Tin were previously published <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Philippakis1" target="_blank">[5]</a>. Motif matches for motifs most relevant to the classification for a given DNA binding domain class: POUHD (V$OCT_01, V$POU1F1_Q6, V$OCT4_02), Tbx (V$TBX5_01, I$BYN_Q6), Myb (V$MYB_Q6), Fkh (V$FOXO3_01, V$FOXO1_Q5, V$FREAC2_01), HD (I$ABDA_Q6, V$CDX5_Q5, V$IFP_03, V$PAX4_02), and Mef2 (V$AMEF2_Q6, V$HMEF2_Q6). These sites were mapped using MAST under default parameters <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Bailey1" target="_blank">[118]</a>. (B) A generic FC enhancer receives differential input from signal-activated, ubiquitous, tissue-restricted and cell type-specific TFs. HD binding motifs are represented as both tissue-restricted and cell type-specific classes since these motifs receive input from both Hox TFs, which are widely expressed in the mesoderm <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Capovilla1" target="_blank">[35]</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Michelson1" target="_blank">[59]</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Enriquez1" target="_blank">[67]</a>, and muscle identity HD TFs—such as Slou, Msh and Ap—which are cell type-specific <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Knirr2" target="_blank">[54]</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Bourgouin1" target="_blank">[68]</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Nose1" target="_blank">[69]</a>. For this diagram, HD binding sites were not subdivided into the distinct binding profiles that have been identified for each individual HD TF (<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Noyes1" target="_blank">[83]</a>, <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1002531#pgen.1002531-Berger1" target="_blank">[126]</a> and B. W. Busser, L. Shokri, S. A. Jaeger, S. S. Gisselbrecht, A. Singhania, M. F. Berger, B. Zhou, M. L. Bulyk and A. M. Michelson, unpublished data).</p

Requirements for NK-2 binding sites for the full activities of multiple tested mesodermal enhancers.

<p>(A) Loss of ßgal reporter (green) in the Lb-expressing SBM (magneta) driven by a version of the <i>lbl</i> enhancer in which the Tin binding sites are inactivated (<i>lbl^noNK-2-lacZ</i>) in stage 14 embryos. Compare to the WT version of the <i>lbl</i> enhancer (<i>lbl^WT-lacZ</i>) in Figure 2A. (B) Normal GFP reporter (green) activity in the <i>ap</i> enhancer in which the Tin binding sites are mutated (<i>ap^noNK-2-GFP</i>) in stage 14 embryos. Compare to the WT version of the <i>ap</i> enhancer (<i>ap^WT-GFP</i>) in Figure 2C. (C) Attenuation of GFP (green) driven by a version of the <i>mib2</i> enhancer in which Tin binding sites are inactivated (<i>mib2^noNK-2-GFP</i>) as compared to ßgal (magneta) driven by a WT version of the <i>mib2</i> enhancer (<i>mib2^WT-lacZ</i>) in stage 12 embryos. Compare to WT versions of both GFP and lacZ reporters in Figure 2E. (D) Loss of GFP (green) driven by a version of the <i>Ndg</i> enhancer in which Tin binding sites are mutated (<i>Ndg^noNK-2-GFP</i>) as compared to ßgal (magneta) driven by a WT version of the <i>Ndg</i> enhancer (<i>Ndg^WT-lacZ</i>) in stage 12 embryos. The ventral <i>Ndg</i> reporter-expressing cells are not in this focal plane but do not express the GFP reporter (data not shown). Compare to WT versions of both GFP and lacZ reporters in Figure 2G.</p