Single cell evaluation of endocardial Hand2 gene regulatory networks reveals HAND2-dependent pathways that impact cardiac morphogenesis

Cardiac morphogenesis is a complex process requiring the synergistic action of multiple tissue types and fine-tuned control of morphogenetic events that are coordinated through the actions of transcription factors within each contributing cell type. Spatial and temporal specific cell signaling between the developing myocardium, the muscular portion of the heart, and endocardium, the inner endothelial lining of the heart, is required for normal cardiogenesis. The embryonic day (E) 12.5 myocardium of the developing ventricles expresses vascular endothelial growth factor A (VEGFA), which signals to the endocardium through its interactions with VEGF receptor 2 (VEGFR2; also known as KDR) to establish the coronary plexus, which will mature to contribute to the coronary arteries (Wu et al., 2012). In the atria, NOTCH signaling initiated within the endocardium communicates with receptors expressed in the myocardium regulating valve and sinoatrial node development (Wang et al., 2020, 2013). Ventricular NOTCH signaling from the endocardium is also required for myocardial BMP10 expression, which is essential for normal trabeculation (Chen et al., 2004; Del Monte-Nieto et al., 2018; Grego-Bessa et al., 2007).

Recent work has revealed that the basic helix loop helix (bHLH) transcription factor HAND2 is required for NOTCH-dependent functions within the endocardium, modulating trabeculation, septation, coronary vascular maturation, as well as endocardial maturation within the embryonic heart (VanDusen et al., 2014a). Conditional endocardial deletion of Hand2 using Nfatc1^Cre (H2CKO) (VanDusen et al., 2014a; Wu et al., 2012) results in embryonic lethality by E12.5. Embryos exhibit tricuspid atresia (TA) or double inlet left ventricle (DILV), in which both tricuspid and mitral valves connect the atria with the left ventricle (LV) (VanDusen et al., 2014a). H2CKOs also exhibit hypoplastic myocardium, an intraventricular septum (IVS) that is shifted to the right generating a smaller right ventricle (RV) and larger LV. H2CKO hearts are hypotrabeculated, and occasionally present with multiple IVS (VanDusen et al., 2014a). In addition to these defects, H2CKO hearts exhibit a pronounced hypervascularization phenotype with an increased number of coronary arteries within the myocardium (VanDusen et al., 2014a).

To investigate transcriptomic changes within H2CKO endocardium, we employed single cell (sc) RNA-seq to identify several gene regulatory networks (GRNs) compromised by HAND2 loss of function. The most notable GRN impacted by the loss of HAND2 is the Apelin Endothelial Signaling Pathway, which is related to shear-stress response. Based on the overlap of significant gene expression changes with established HAND2 DNA occupancy data (Laurent et al., 2017), we sought to identify putative HAND2-dependent endocardial transcriptional enhancers. We selected five genes that showed robust expression changes associated with bound HAND2 to further evaluate: Igf2, Igf2R, Ptn, Tmem108 and Klf2. Putative enhancer sequences for two of the selected genes, Igf2R and the shear-stress master regulator gene Klf2, exhibited functional endocardial/endothelial transcriptional enhancers via F0 transgenic reporter analysis. We further characterized the activity of a 1.8 kb Klf2 enhancer located −50 kb upstream of the Klf2 transcriptional start site (TSS). We observed that the −50 kb Klf2 enhancer is active within the early developing vasculature endothelium and, importantly, within the endocardium, recapitulating the endogenous Klf2 expression pattern. To determine whether HAND2 is necessary for Klf2 enhancer endocardial expression in vivo, we interrogated the activity of the −50 kb Klf2 enhancer on the H2CKO background. Indeed, our data revealed that, in the absence of HAND2, activity of the −50 kb Klf2 enhancer is robustly reduced within trabecular endocardium but is unaffected within the systemic vasculature and developing valves. Gene edited deletion of the −50 kb Klf2 [Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb)}] enhancer demonstrated its requirement for Klf2 endocardial expression. Collectively, these findings demonstrate several previously unreported endocardial HAND2-dependent gene regulatory pathways, including the shear-stress response pathway, mediated in part through HAND2 regulation of Klf2.

To further investigate the role of HAND2 within the developing endocardium, we crossed Hand2-conditional mice (H2^fx/fx;R26R^mTmG/mTmG) (Morikawa et al., 2007; Muzumdar et al., 2007) with the endocardial-specific Nfatc1^Cre (Wu et al., 2012) to generate H2CKO (Nfatc1^CreHand2^fx/fxR26R^mTmG/wt) as well as control littermates (Hand2^fx/+R26R^mTmG/wt) and isolated E11.5 hearts for scRNA-seq analysis using the 10x Genomics platform. In H2CKO hearts, Cre-recombinase mediated recombination led to switching of tdtomato epifluorescence to GFP epifluorescence, which allowed for the quick identification of Cre-positive embryos. Rapid PCR genotyping was used to identify the Hand2-conditional allele status of the embryos. We sequenced 13,885 unique barcodes from a single H2CKO heart, and 14,259 barcodes from a single control heart. Based on the high expression of hemoglobin genes, we excluded 5828 barcodes from H2CKO and 3150 barcodes from control hearts. Next, we excluded barcodes where the total number of genes was greater than 2500 (indicating multiplets) from the analysis. The remaining 6232 barcodes from H2CKO and 5408 barcodes from control hearts were used for further analysis. Non-linear dimensionality reduction using uniform manifold approximation and projection (UMAP) plots resulted in 13 transcriptionally distinct clusters (Fig. 1A; Table S1). Fig. S1 displays control and H2CKO cells mapped separately.

Cluster identity was assigned by comparing gene expression of a gene in the control cluster against expression of the same gene in all other control clusters combined at a threshold of 0.25 log2FC, which establishes a rigorous threshold for significance (Table S1). Clusters 0 (red) and 1 (pink) represented 2793 cells, exhibited similar transcriptional profiles and expressed the following cardiomyocyte marker gene transcripts: Myh6 (99.4% cells in cluster 0, 98.1% cells in cluster 1), coding for alpha myosin heavy chain (αMHC), and actinin alpha 2 (Actn2; 99.9% of cells in clusters 0 and 1; Fig. 1A; Fig. S2; Table S1). The 1030 cells within cluster 2 (brown) expressed the extracellular matrix protein coding gene periostin (Postn, 99.8% of cells), the bHLH transcription factor gene Twist1 (99.1% of cells), the matricellular protein coding gene transforming growth factor beta induced (Tgfbi, 98% of cells) and represent atrioventricular (AV) cushion cells that have undergone endothelial-to-mesenchymal transition (EMT, Fig. S2; Table S1). Cluster 3 (orange) consisted of 993 cells that were identified as outflow tract mesenchyme based on the expression of the neurovascular guiding factor semaphorin 3c (87% of cells) and bone morphogenetic protein 4 (Bmp4, 88% of cells; Fig. S2; Table S1). Clusters 4 (light purple; 740 cells) and 5 (dark purple; 686 cells) consisted of cardiac neural crest cells as determined by the expression of Twist1 (99.8% of cells in cluster 4, 100% in cluster 5), insulin growth factor (Igf1, 98% of cells in both clusters), and high mobility group transcription factor Sox9 (94% of cells in cluster 4 and 96% of cells in cluster 5, Fig. S2; Table S1). The 594 cells in cluster 6 consisted of epicardial cells (light green), undergoing transition to a fibroblast phenotype as marked by the expression of the growth factor pleiotrophin (Ptn, 100% of cells), T-box transcription factor Tbx18 (90% of cells) and bHLH transcription factor Tcf21 (98% of cells; Fig. S2; Table S1). Cluster 7 consisted of 564 cells (light blue) that were identified as endocardial cells by the expression of the transmembrane transport protein Ramp2 (100% of cells), the vascular endothelial cadherin (Cdh5, 99% of cells, Fig. S2; Table S1), platelet endothelial cell adhesion molecule (Pecam1, 95% of cells), and the nuclear factor of activated T cells (Nfatc1, 64% of cells). Cluster 8 consisted of 542 cells and was identified as a second epicardial cell cluster (dark green) expressing Ptn (100% of cells), Tbx18 (97% of cells) and Tcf21 (97% of cells; Fig. S2; Table S1). Cluster 9 consisted of 462 cells and represented a second endocardial cell cluster (dark blue) expressing Ramp2 (100% of cells), Cdh5 (95% of cells), Pecam1 (98% of cells) and Nfatc1 (88% of cells; Fig. S2; Table S1). Cluster 10 did not express any gene that exhibited a log2 fold change of expression within the cluster and outside the cluster greater than 0.25, thus this cluster (grey) has remained undefined. Cluster 11 represented the conduction system cell cluster (yellow) marked by expression of calcium channel, voltage-dependent, α2/δ2 subunit 2 (Cacna2d2, 97% of cells) and calcium channel, voltage-dependent, T type, α 1H subunit, (Cacna1h, 65% of cells). Cluster 12 (black) represented the lymphocyte population as indicated by expression of interferon induced transmembrane protein 3 (Ifitm3, 47% of cells) and histocompatibility 2, D region locus 1 (H2-D1, 93% of cells; Fig. S2; Table S1). Nfatc1 expression robustly marked both cluster 7 (64% of cells) and cluster 9 (88% of cells; Fig. 1B). Comparison of H2CKO and control barcodes specific for Hand2 expression exhibited robust downregulation within these two endocardial clusters in the presence of Nfatc1^Cre (cluster 7 log2FC −1.5, P=1.5×10⁻⁵²; cluster 9 log2FC −1.66, P=9.3×10⁻⁵⁵; Fig. 1B; Fig. S3). Note that cluster 10 (undefined) maintained Hand2 expression post-deletion (Fig. 1B; Fig. S3).

To examine transcriptome data in an unbiased fashion, we employed Ingenuity Pathway Analysis (IPA) on differentially expressed genes (log2FC>±0.5) from H2CKO and control cells (Table 1; Fig. S4; Table S2). Loss of Hand2 within the endocardium led to significant changes in developmental, morphological, and cardiovascular gene regulatory networks (Table 1; Fig. S5). IPA on cluster 7 indicated that Hand2 is downregulated within the Cardiac Hypertrophy Signaling (Enhanced) canonical pathway (z-score −1.807, −log P-value 2.8; Table S2), which is close to a significant z-score absolute value of 2.

We also observed abnormal expression of a number of endocardial transcripts including endothelin converting enzyme 1 (Ece1), a shear-stress responsive gene expressed by vascular endothelial cells and required for formation of patent vasculature in the developing heart (Masatsugu et al., 2003; Robinson et al., 2014). Ece1 was significantly downregulated within H2CKO endocardium (cluster 7 log2FC −1.16, P=2.69×10⁻²⁸; cluster 9 log2FC −0.45, P=9.05×10⁴; Table S3). Concomitantly, the Ece1 substrate endothelin 1 (Edn1), a potent vasoconstrictor that is secreted by endothelial cells when laminar flow induces shear stress (Morawietz et al., 2000), was upregulated within H2CKO endocardium (cluster 7 log2FC 0.58, P=1.9×10⁻⁹). Previous work has shown that EDN1 signaling lies upstream of Hand2 within the cranial neural crest cells during craniofacial morphogenesis (Charité et al., 2001; Clouthier et al., 2000). Thus, the observed increase in Edn1 could reflect a feedback compensation as the result of the endocardial loss of Hand2.

Fibronectin 1 (Fn1), a component of the extracellular matrix secreted by endothelial cells, was significantly downregulated within clusters 7 and 9 (cluster 7 log2FC −1.17, P=1.4×10⁻²¹; cluster 9 log2FC −0.66, P=1.6×10⁻⁴, Table S3). Fn1 appeared in multiple IPA pathways (Table 1) including: Wound Healing Signaling Pathway (z-score −3.606, −log P-value 4.67), Pulmonary Fibrosis Idiopathic Signaling Pathway (z-score −3.273, −log P-value 9.66) and Tumor Microenvironment Pathway (z-score −2.449, −log P-value 2.1).

The mechanosensitive transcription factor hypoxia inducible factor 1 alpha (Hif1a) was significantly downregulated within endocardial clusters (cluster 7 log2FC −0.58, P=4.5×10⁻¹²; cluster 9 log2FC −0.48, P=2.3×10⁻¹⁰; Table S3; Feng et al., 2017). IPA analysis revealed that related canonical pathways which included Hif1a were disrupted (Table 1; Table S2): Pulmonary Healing Signaling Pathway (z-score −2.65, −log P-value 1.87) and Tumor Microenvironment Pathway (z-score of −2.449, −log P-value 2.1). The Tumor Microenvironment Pathway also included Igf2, the gene that encodes insulin-like growth factor 2, which was significantly downregulated in H2CKO endocardial clusters (cluster 7 log2FC −1.44, P=3.2×10⁻⁴⁶; cluster 9 log2FC −1.81, P=2.1×10⁻⁶⁵). Indeed, Hif1a transcriptionally regulates Igf2 via hypoxia responsive elements at the Igf2 locus (Feldser et al., 1999). Within endocardial cells, angiogenesis requires the action of the shear-stress master regulator KLF2 (Nigro et al., 2011).

Given the presence of two endothelial/endocardial clusters and to better understand the differences within these cells, we undertook the direct comparison of gene expression within clusters 7 and 9 between the control and H2CKO cell populations (Fig. S6; Table S4). Results showed that 785 genes were differentially expressed between control in clusters 7 and 9. These genes included Klf2, (log2FC −0.536095853), Hey2 (Seya et al., 2021) (log2FC −0.610741244), the NOTCH-dependent ligand Wnt4 (Luxán et al., 2016) (log2FC −0.857926425) and the endocardial specific Irx6 (Mummenhoff et al., 2001) (log2FC 0.288606862; Table S4, tab control). Similar analysis on H2CKO data revealed that 981 genes were differentially expressed between clusters 7 and 9, including the aforementioned examples (Table S4, tab H2CKO). In addition, 501 common genes were differentially regulated between control and H2CKO clusters 7 and 9, with 284 genes uniquely regulated in controls and 480 genes uniquely regulated in H2CKO (Fig. S6; Table S4, tabs control only and H2CKO only). Given that there was not a significant amount of coronary vasculature present at E11.5 (Ivins et al., 2015), it is possible that these two clusters represented distinct populations of maturing endocardium, where cluster 7 cells might contribute to the future coronary vasculature as H2CKO hearts exhibited hypovascularized ventricles (VanDusen et al., 2014a).

IPA also revealed that H2CKO mutants exhibited a significant downregulation of the Apelin Signaling Pathway, with a z-score of −2 (−log P-value 1.47, Table 1). Apelin is an angiogenic factor that controls migration of endothelial cells and is required for the normal development of blood vessels (Helker et al., 2020; Kwon et al., 2016; Lu et al., 2017). The Apelin Signaling IPA pathway includes a major contributing factor to normal vasculogenesis and ventricular morphogenesis within the embryonic heart – shear-stress signaling (Haack and Abdelilah-Seyfried, 2016). We observed that the gene coding for the shear-stress regulated transmembrane receptor heart of glass (Heg1) was significantly downregulated within both endocardial clusters (cluster 7 log2FC −0.36, P=6.35×10⁻⁹; cluster 9 log2FC −0.9 P=1×10⁻¹⁵; Table S3). Heg1 zebrafish mutants exhibit significant vascular malformations (Kleaveland et al., 2009). Interestingly, the transcription factor KLF2 is a direct regulator of Heg1 expression (Razani et al., 2001; Zhou et al., 2016). KLF2 is a well-studied shear-stress response transcription factor considered the master regulator of this response (Bhattacharya et al., 2005; Chiplunkar et al., 2013; Sangwung et al., 2017).

As the hypervascularization phenotype observed in the H2CKO mutants could be caused by defective angiogenesis, and KLF2 is a major regulator of angiogenesis, we employed in-situ hybridization to closely examine both Hand2 and Klf2 transcripts within E12.5 control and H2CKO embryo hearts (Fig. 2). Results showed that Hand2 expression within the trabecular endocardium of the ventricles was significantly reduced in H2CKO hearts when compared with controls (Fig. 2A,C). Klf2 was also robustly expressed within the ventricular endocardium (Fig. 2B,B′, black arrowheads), and particularly within areas of high shear stress such as the endocardial lining of the AV canal (Fig. 2B′, red arrowheads; Chiplunkar et al., 2013). We observed that Klf2 expression within the ventricular endocardium of H2CKO hearts was greatly reduced; however, Klf2 expression within the endocardial lining of the AV canal as well as within the systemic vasculature was maintained (Fig. 2D,D′, compare tissues marked by red and black arrows). Next, we examined genes downstream of KLF2 that were significantly changed within H2CKO endocardial clusters 7 and 9 (Table S5). Indeed, we observed significant changes in gene expression within KLF2 target genes, which suggested that loss of HAND2 in the endocardium reduced Klf2 expression as well the expression of Klf2 target genes (Table S5).

The differential gene expression profiles in endocardial clusters from wild-type (WT) and H2CKO hearts suggested a direct interaction of the HAND2 transcription factor with cis-regulatory elements in the respective loci for transcriptional control. We selected a subset of genes (Igf2, Tmem108, Ptn, Igf2R and Klf2) for which distinct HAND2 interaction peaks were identified in the respective regulatory domains by determining regions of evolutionary conservation and HAND2 DNA binding (Fig. 3; Table 2). For identification of HAND2 target regions, we used established E10.5 chromatin immunoprecipitation (ChIP)-seq data from mouse embryonic hearts expressing a Hand2^3xFlag knock-in allele (Laurent et al., 2017). To validate and define cardiac in vivo activities of putative enhancer regions in the selected subset of gene loci we employed lacZ transgenic reporter assays, involving enSERT, a method for CRISPR-mediated site-directed reporter transgenesis targeting the H11 locus (Kvon et al., 2020; Osterwalder et al., 2022).

Three of our putative HAND2-dependent enhancers that exhibit evolutionary conservation and HAND2 DNA occupancy did not exhibit endocardial/endothelial enhancer activity (Fig. S7). IGF2 is a secreted growth factor expressed within the epicardium and endocardium during heart development (Shen et al., 2015) and is highly downregulated within clusters 7 and 9 in H2CKO mutants (Table S3). At the Igf2 genomic locus, a conserved non-coding element (CNE) located 70 kb 3′ to the coding region (Fig. S7A) showed pronounced HAND2 DNA occupancy. However, analysis by enSERT exhibited no enhancer activity of this element (n=6/9 tandems; Fig. S7B,B′). TMEM108 is implicated as a marker for progenitor epicardial cell populations although its role within the endocardium is currently unclear (Bochmann et al., 2010). Tmem108 is significantly downregulated within endocardial clusters 7 and 9 (Table S3); however, a HAND2-occupied CNE located 5′ of Tmem108 exhibited no enhancer activity (n=2/3 tandems; Fig. S7C-D′). Ptn codes for a secreted cytokine, is an inducer of EMT, and is mitogenic to endothelial cells, resulting in angiogenesis (Perez-Pinera et al., 2008). Ptn is significantly downregulated within both clusters 7 and 9 (Table S3). We tested a HAND2-occupied CNE located 3′ of Ptn; however, results showed no E11.5 heart expression (n=2/3 tandems; Fig. S7E-F′).

Three of our putative HAND2-dependent enhancers that exhibited robust HAND2 DNA occupancy also exhibited endocardial/endothelial enhancer activity. We successfully interrogated a HAND2-occupied CNE located 21 kb 5′ of the Igf2R TSS (Fig. 3A). Igf2R is robustly expressed within the endocardium (Wang et al., 2019) and lacks a tyrosine kinase domain acting as a negative regulator of IGF2 (Braulke, 1999; Ludwig et al., 1996). Igf2r was significantly downregulated within H2CKO endocardial clusters 7 and 9 (Table S3). EnSERT analysis of the Igf2R HAND2-occupied CNE resulted in 7/8 transgenic embryos (n=5 tandems) with cardiac-specific lacZ staining at E11.5, including the endocardium, thus uncovering a novel endocardial enhancer (Fig. 3B,B′).

From the observation that both Igf2R and Igf2 are downregulated within the endocardial clusters, we wanted to determine whether endocardial proliferation was affected in H2CKO mutants. We conducted differential abundance analysis to determine cell type representation within H2CKO and control barcodes (Table S6; Fig. S8). Results indicated a significant increase in the number of barcodes in the endocardial population (cluster 7), as well as an increased number of barcodes in the cardiomyocyte population (clusters 0 and 1), which suggested that the loss of IGF2R might have affected the cell numbers in H2CKO mutants. The discovery of a HAND2-binding CNE −21 kb 5′ of the Igf2r TSS that drove endocardial-specific reporter expression supported this idea (Fig. 3A-B′). In order to determine the evolutionary conservation of this element, we used CLUSTWAL analysis, which indicated that conservation was limited within mammals (Fig. S9).

We next interrogated HAND2-occupied putative cardiac enhancers within the locus of the shear-stress master regulator KLF2. Expression analysis in both the scRNA-seq analysis and in situ hybridization (ISH) experiments revealed dynamic Klf2 regulation within the heart (Fig. 2B,B′,D,D′). Klf2 expression was specifically downregulated by HAND2 within the ventricular endocardium; however, Klf2 endothelial expression was maintained within AV cushion endocardium (Fig. 2B′,D′). This observation was consistent with the downregulation of Hand2 expression within the AV cushion endocardium post EMT, as previously reported (VanDusen et al., 2014a), and suggested that there was also HAND2-independent regulation of Klf2 transcription within the AV cushion endocardium, as well as the systemic vasculature, in which Hand2 is not expressed (VanDusen et al., 2014a).

Previously, a CNE located 100 bp upstream of the Klf2 TSS had been shown to be responsive to shear stress (Huddleson et al., 2004), but HAND2 DNA occupancy data did not indicate HAND2 DNA-binding within this element (Laurent et al., 2017). Interestingly, we identified two HAND2-occupied CNE within Klf2, one located −16 kb upstream of the Klf2 TSS and another located −50 kb relative to the Klf2 TSS (Fig. 3C-F′). The −16 kb Klf2 CNE exhibited a high level of HAND2 occupancy and sequence conservation (Fig. 3C). Results showed that two out of eight F0 transgenics at E11.5 exhibited β-galactosidase staining within embryonic vasculature and endocardium (Fig. 3D,D′). Out of these two, we observed only one embryo that showed consistent staining throughout the left and right ventricular endocardium. The −50 kb Klf2 CNE also exhibited robust HAND2 DNA occupancy and sequence conservation (Fig. 3E), and therefore was used to generate E11.5 F0 transgenics as well. Out of seven transgenics generated, four exhibited robust β-galactosidase staining within the endocardium and systemic vasculature (Fig. 3F,F′).

Given the higher consistency in endocardial/endothelial activity that we observed from the −50 kb Klf2 enhancer element (57% of F0s), we generated stable lacZ transgenic lines using the −50 kb CNE (Fig. 4A). At E11.5, out of the 12 transgenic lines generated, five (41%) recapitulated consistent and robust β-galactosidase staining within the endocardium and systemic vasculature (Fig. S10). The remaining seven lines exhibited no observable β-galactosidase staining at E11.5. We next examined additional embryonic time points using a single line (line #901, Fig. S10). Analysis of reporter activity in E7.5 embryos revealed β-galactosidase staining within endothelial precursors, the blood islands (Fig. 4B), and within the dorsal aorta at E8.5 (Fig. 4C). At E9.5 and E10.5, the −50 kb CNE robustly drove β-galactosidase expression within endothelial structures, within the branchial arches and within intersomitic blood vessels (Fig. 4D,E arrow). Histological cross sections of β-galactosidase-stained torsos counterstained with nuclear fast red (NFR) at E10.5 revealed endocardial-specific staining at this time point (Fig. 4F-G′). Thus, we concluded that the −50 kb Klf2 CNE functions as a transcriptional enhancer within the endothelial cells of the developing embryonic vasculature, endocardium and AV cushions (Fig. 4F′,G′).

Motif analysis of the −50 kb Klf2 enhancer revealed the presence of three conserved Eboxes (Fig. 4H; Fig. S11) that lie within this established HAND2 DNA occupancy peak (Laurent et al., 2017). In order to determine whether HAND2 was able to interact with any of the three conserved Eboxes within this Klf2 enhancer, we conducted ChIP assays in NIH-3T3 cells by co-transfecting plasmids encoding a 5′ Myc-tagged Hand2 and an untagged E12 (Tcf3) (Fig. 4H′). Negative controls used pCS2+myc samples immunoprecipitated with and without αMyc, and Myc-Hand2 immunoprecipitated without αMyc. Using ChIP-PCR, we were able to observe HAND2 DNA binding at the 5′-most Ebox (Ebox1 CACCTG) within the Klf2 enhancer in a dose-dependent manner (Fig. 4H′). The controls for this experiment employed primers recognizing the mouse Rpl30 gene (Fig. 4H″). Taken together, the in vitro data supported HAND2 directly binding to and transcriptionally regulating Klf2 through the −50 kb Klf2 CNE.

To assess whether the endocardial −50 kb Klf2 enhancer was dependent on HAND2 in vivo, we crossed the −50 kb Klf2 enhancer lacZ transgenic reporter with the endocardial-specific H2CKO. E11.5 embryos were β-galactosidase stained, sectioned and counterstained with NFR (Fig. 4I). The −50 kb Klf2 enhancer embryos (lacZ+ Hand2^fx/+) exhibited positive staining of trabecular endocardium (Fig. 4I,I′). In comparison, −50 kb Klf2 enhancer lacZ reporter H2CKO embryos (lacZ+; Nfatc1^cre Hand2^fx/fx) showed a robust reduction in ventricular endocardial β-galactosidase staining (Fig. 4J,J′, arrows), whereas β-galactosidase staining within the endocardium over the developing AV cushions and systemic vasculature was maintained (Fig. 4J).

The −50 kb Klf2 enhancer activity recapitulated the Klf2 mRNA expression pattern throughout the embryo, with activity within the developing vasculature and the endocardium; however, crossing the −50 kb Klf2 enhancer lacZ reporter to the H2CKO background only altered enhancer β-galactosidase staining within the ventricular endocardium and did not appreciably alter expression within the systemic vasculature or within the AV cushions. We found this result to be completely consistent with ISH analysis of Klf2 expression in the H2CKO mutants (Fig. 2B′,D′). These data suggested to us that HAND2 DNA binding within this −50 kb Klf2 endothelial/endocardial enhancer was necessary for its activity within the ventricular endocardium.

As the −50 kb Klf2 endothelial/endocardial enhancer recapitulated Klf2 vascular expression, we next tested its requirement for maintenance of Klf2 expression via CRISPR-mediated genomic deletion in mice (Fig. S12A). Eleven enhancer-deleted lines were obtained and four of these were crossed for two generations with WT mice before they were intercrossed for homozygosity. Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb)} mice were viable and were born at Mendelian frequencies in four outcrossed founder lines (Fig. S12B). A single line was then set up for timed pregnancies and E11.5 embryos were evaluated for Klf2 expression (Fig. 5A,B). Klf2 expression was visibly lower within the endocardium. In contrast, Hand2 expression within adjacent sections appeared to be unchanged between controls and Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb)} hearts (Fig. 5C,D). We performed qRT-PCR on E11.5 ventricles to confirm that the observed Klf2 expression drop in the Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb)} homozygous hearts was significant. We isolated eight Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb)} and ten WT ventricles for qRT-PCR (Fig. 5E). As predicted by the ISH analysis, Klf2 expression levels were significantly lower (P<0.001) in Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb} ventricles compared with control ventricles, ∼60% of what is observed in WT (Fig. 5E). Results showed that Hand2 expression was unchanged within Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb)} hearts compared with WT controls (Fig. 5E). Given that the −50 kb enhancer recapitulated all Klf2 embryonic expression but was only affected by HAND2 within the ventricular endocardium, the observed 40% decrease in endocardial expression was in line with our observations.

Loss of Hand2 within the endocardium disrupts NOTCH signaling resulting in a hypotrabeculated single ventricle composed of hypervascularized free walls (VanDusen et al., 2014a). To gain a better understanding on the gene regulatory networks in which HAND2 facilitates ventricular morphogenesis downstream of NOTCH1, we used scRNA-seq at E11.5, combined with established HAND2 DNA occupancy data (Laurent et al., 2017) to interrogate the role of HAND2 in regulating the endocardial gene regulatory networks. IPA analysis revealed a number of pathways known to be required for heart development that showed misregulation within the identified endothelial/endocardial cell populations (Fig. 1). These analyses revealed disruption in several pathways such as wound healing, pulmonary fibrosis and healing, tumor microenvironment (including HIF1α signaling) as well as the Apelin pathway, which includes shear-stress response regulation and is the pathway most relevant to endocardial roles in cardiogenesis (Table 1; Fig. S5; Table S2). Collectively, these pathways play roles in the endocardial response to vascularization of the myocardium, organ growth and communication with the underlying myocardium coordinating septation and trabeculation. A number of significantly regulated genes, for example Fn1, Ece1 and Edn1, exhibit altered expression, but do not exhibit robust HAND2 DNA occupancy in cis (Laurent et al., 2017). Although such genes are influenced by HAND2 function, they are likely not transcriptionally regulated by HAND2 directly; nevertheless, their altered expression fits with HAND2 function in previous studies. In epicardial Hand2 deletion, although Fn1 expression is unaltered when comparing control with mutants, FN1 organization is altered within H2CKO epicardial cells (Barnes et al., 2011). During jaw morphogenesis, Hand2 has been established as lying downstream of EDN1 signaling and plays an important negative feedback role once activated, by repressing Dlx5 and Dlx6 expression within the ventral-most portion of the mandible mesoderm (Barron et al., 2011; Charité et al., 2001; Clouthier et al., 2000; Vincentz et al., 2016).

We chose five target genes – Igf2, Igf2R, Ptn, Tmem108 and Klf2 – to investigate further for putative endocardial/endothelial HAND2-dependent enhancers (Fig. 3; Fig. S7), based on our comparisons of highly misregulated genes with robust HAND2 DNA occupancy data to locate potential cis-regulatory elements (Laurent et al., 2017). CNE peaks bound by HAND2 from Igf2, Tmem108 and Ptn did not reveal any transcriptional activity (Fig. S7). More interestingly, we discovered three endocardial/endothelial enhancers that did have transcriptional activity, one CNE 5′ of the Igf2R TSS, and two CNE upstream of the Klf2 TSS (Fig. 3).

The crucial source of IGF2 in the heart is from the epicardium (Shen et al., 2015). Epicardial IGF2 diffuses into the heart where it can bind to its receptors, including IGF2R. Binding to IGF2R facilitates IGF2 degradation within the lysosomes (Harris and Westwood, 2012). Knockout of IGF2R within endothelial cells using Tie2Cre does not result in embryonic lethality; however, cardiac-specific phenotypes have not been examined (Sandovici et al., 2022). Our data suggest that IGF2R plays a role within cardiac endothelium in a HAND2-dependent manner. We observed a significant change in cell numbers within H2CKO mutants and controls (Table S5), although it is yet to be determined whether increased proliferation is the cause. Cell proliferation within E10.5 right ventricle of Tie2Cre-mediated H2CKO hearts did not show a significant difference in cell numbers (VanDusen et al., 2014a,b).

As KLF2 is a known master regulator of shear-stress response and is a significant regulator within the Apelin regulatory network, we engineered a stable Klf2 reporter line using the more robustly consistent endothelial/endocardial CNE located at −50 kb of the Klf2 TSS. Reporter expression analysis reveals that this Klf2 CNE recapitulates all the Klf2 endothelial/endocardial expression and is dependent on HAND2 only within the endocardium, correlating directly with our Klf2 mRNA expression regulation data (Figs 2 and 4). KLF2 is expressed in regions of endothelium exposed to high shear stress (Goddard et al., 2017). In the developing heart such high shear-stress regions include the endocardium overlying the developing valves and the developing ventricular trabeculae, with Klf2 expression levels varying within regions of differential shear force (Goddard et al., 2017). Loss of HAND2 does not appear to impact Klf2 expression within the regions of endocardium overlying the developing valve cushions where Hand2 expression is already downregulated (VanDusen et al., 2014a,b). Previous work characterizing conserved non-coding elements at the Klf2 genomic locus identified a 60 bp enhancer element located 100 bp upstream of Klf2 TSS that is responsive to shear stress within mouse endothelial cells in culture (Huddleson et al., 2004). It is currently unclear whether either the −16 or −50 kb CNE enhancers are shear-stress responsive, but it is clear that the −50 kb enhancer can recapitulate all Klf2 mRNA expression domains during embryogenesis and the majority of its activity is HAND2-independent, given that HAND2 is not expressed within the systemic vasculature (Fig. 4). As one would expect, the −50 kb Klf2 enhancer element contains other consensus binding sequences including the myocyte enhancer factor 2 (MEF2) family of transcription factors that are established regulators for vascular homeostasis and are transcriptional activators of Klf2 (De Val and Black, 2009; Lu et al., 2021). Analysis of DNA occupancy data shows that both the −16 kb and −50 kb Klf2 enhancers have conserved MEF2C binding (Akerberg et al., 2019). Given the established role of MEF2C in endothelial integrity and homeostasis, it is not surprising that the loss of Hand2 does not lead to loss of vascular KLF2 expression.

Lineage tracing analysis shows that the endocardium is a primary source of cells that eventually gives rise to coronary vessels (Sharma et al., 2017). Studies in mouse models demonstrate that both endocardial and epicardial cells migrate into the myocardium to give rise to patent vessels (Sharma et al., 2017) and that these coronaries form within different zones of the myocardium (septum versus free walls of the ventricles) (Chen et al., 2014). This suggests that coronary angiogenesis is driven by distinct mechanisms within different regions of the developing heart. The cellular origin of coronary vasculature is a source of some debate, the current consensus being that coronaries of the ventricular free wall are derived from the epicardium and the sinus venosus, whereas interventricular septal coronaries are derived from the ventricular endocardium (Phansalkar et al., 2021; Rhee et al., 2021; Zhang et al., 2016).

Klf2 undergoes a robust shear-stress response, as at least 50% of the highly regulated flow genes are dependent on the upregulation of Klf2 (Parmar et al., 2006). Klf2 expression within the endocardial cells of the ventricular wall fated to contribute to coronaries in these endocardial cells is HAND2 dependent. Indeed, one of the most striking observations in H2CKO heart endocardium and vasculature is the persistent expression of Lyve1 beyond its normal endocardial downregulation by E13.5 (VanDusen et al., 2014a). LYVE1-expressing endocardium ultimately contributes to peripheral cardiac macrophages and the developing lymphatic vasculature of the heart, in which vessel pressures are far less than those encountered in blood vasculature (Pinto et al., 2012). It is an appealing idea that a defective shear-stress response of the ventricular endocardium could result in improper development and/or maturation of the ventricular endocardium into the correct sub-fates that result in hypervascularization of the ventricular walls composed of an immature, more lymphatic-like, endothelium. Further support for this idea comes from multiple lines of evidence demonstrating that KLF2 inhibits angiogenesis by interacting with the Kdr promoter (Bhattacharya et al., 2005), as the loss of KLF2 also leads to hypervascularization (Kawanami et al., 2009). In our H2CKO data, we observed a modest increase in Kdr expression in cluster 7 (log2FC 0.13, not significant) which could be supportive of this possible mechanism.

Multiple genetic knockouts have been generated to study KLF2 function within endothelial cells. The Klf2 systemic knockout is embryonically lethal between E12.5 to E14.5 due to severe intra-embryonic and intra-amniotic hemorrhaging (Kuo et al., 1997; Wani et al., 1998). Endothelial knockout of Klf2 (and the related Klf4) using tamoxifen-inducible Cdh5-Ert2Cre in 8- to 10-week-old adult mice causes vascular leakage leading to hemorrhaging and death (Sangwung et al., 2017). Embryonic endothelial knockout of Klf2 using Tie2Cre exhibits increased systolic stroke volumes and high output heart failure leading to death at E14.5 due to abnormal vessel tone (Lee et al., 2006) and endocardial knockout of Klf2 using Nfatc1^Cre results in embryonic lethality by E14.5 due to septal defects arising from the failure of cushion remodeling (Goddard et al., 2017). Moreover, work in zebrafish demonstrates that flow-responsive Klf2 activates notch signaling, through a mechanism that employs endocardial primary cilia (Li et al., 2020).

Given that Klf2^{Δ-50:(3.9kb)/Δ-50:(3.9kb)} mice exhibit only a 40% reduction in Klf2 endocardial expression and appear to maintain systemic vascular expression through other identified enhancers (Fig. 5), it is not surprising that the removal of this −50 kb Klf2 CNE does not result in embryonic lethality and that mice are viable and fertile. What we do not know currently is the Klf2 expression threshold that results in the observed embryonic vascular phenotypes or whether any KLF2 endocardial-specific phenotypes contribute to the observed embryonic lethality. Collectively, these data demonstrate that HAND2 integrates endocardial transcriptional networks reaching beyond the NOTCH pathway and including shear-stress response, thereby revealing a number of important roles during endocardial morphogenesis.

Hand2^fx/fx mice (Morikawa and Cserjesi, 2008; The Jackson Laboratory, strain 027727) and Nfatc1^cre (Wu et al., 2012) were genotyped as described previously (VanDusen et al., 2014a). The University of Michigan Transgenic Animal Model Core generated lacZ transgenic enhancer lines in the FVB background. Twelve transmitting founder lines were screened for X-gal staining and enhancer activity. Transgenic founders and embryos were genotyped using primers spanning the enhancer and HSP68 promoter 5′-AGCCTGTGAGAGAGACCCAT-3′ and 5′-GATGTTCCTGGAGCTCGGTA-3′. Genotyping for other alleles was carried out using Southern blots as previously described (George and Firulli, 2021). All animal maintenance and procedures were performed in accordance with the Indiana University School of Medicine protocol 20090, and University of Michigan School of Medicine. Animal work at Lawrence Berkeley National Laboratory (LBNL) was reviewed and approved by the LBNL Animal Welfare Committee.

E11.5 embryos were dissected in cold PBS and placed in PBS with 1% fetal bovine serum (FBS) solution on ice until dissociation (∼3 h). Yolk-sac DNA was extracted (QuickExtract DNA Extraction Solution, Epicentre) and used for genotyping to distinguish heterozygous and homozygous Hand2 conditional alleles. The Rosa^mTmG allele fluorescence was used to determine Nfatc1^cre status. Dissected cardiac tissue was incubated in 750 μl TrypLE (Thermo Fisher Scientific) for 5 min at 37°C, triturated with a 200-μl wide-bore pipette tip. The cell suspension was quenched with 750 μl DMEM with 10% FBS. Cells were filtered through a 30-μm cell strainer (MACS SmartStrainer), centrifuged at 300 g for 5 min, and washed once with 750 μl PBS with 0.5% bovine serum albumin (BSA). Cells were resuspended in 30 μl PBS with 0.5% BSA (10x Genomics). Single-cell droplet libraries from this suspension were generated using the Chromium NextGEM Single Cell 3′ Reagent Kits User Guide, CG000204 Rev D (10x Genomics), according to the manufacturer's instructions. Briefly, each clean single cell suspension was counted with hemocytometer under a microscope for cell number and cell viability. Only single cell suspensions with a viability of >90% and minimal cell debris and aggregation were used for further processing. The resulting library was sequenced in a custom program for 28b plus 91b paired-end sequencing on Illumina NovaSeq 6000. About 50,000 reads per cell were generated and 91% of the sequencing reads reached Q30 (99.9% base call accuracy).

Sequenced reads were aligned to a mouse transcriptome reference built from GRCm38.p6 (Genome Reference Consortium Mouse Build 38 patch release 6) combined with eGFP and dTomato gene sequences using the software 10x Genomics Cell Ranger 5.0.1 (Zheng et al., 2017). Reads from the cells associated with more than 1000 unique molecular identifiers from hemoglobin-related genes (Hbb-bt, Hbb-bs, Hbb-bh2, Hbb-bh1, Hbb-y, Hba-x, Hba-a1, Hbq1b, Hba-a2 and Hbq1a) were excluded from further analysis. The downstream data exploration and differential gene expression analysis was conducted using the R package, Seurat V4 (Hao et al., 2021). As per the standard pre-processing workflow for scRNA-seq data in Seurat, cells with more than 2500 unique features were filtered out. The feature expression values for each cell were normalized using the standard ‘LogNormalize’ method with default parameter values. The Seurat objects derived from control (WT) and H2CKO mutant data were integrated using the anchors found using canonical correlation analysis (CCA) with the neighbor search space specified using 1 to 20 dimensions [FindIntegrationAnchors(reduction=“cca”, dims=1:20)]. The integrated dataset was subjected to linear transformation followed by linear dimensionality reduction using principal component analysis (PCA). Clusters were identified from the shared nearest neighbor graph [FindNeighbors(reduction=“pca”, dims=1:20)] with the resolution set to 0.5 [FindClusters(resolution=0.5)] and were visualized using the UMAP non-linear dimensional reduction technique. For each cluster, the differentially expressed genes between control and H2CKO genotypes were called using a Wilcoxon Rank Sum test [FindMarkers(test.use=“wilcox”)]. Genes with Bonferroni corrected P-values not more than 0.05 were considered significantly differentially expressed. For IPA analysis, the pathways relevant to the significantly differentially expressed genes (FDR≤0.05) were identified using the Core Analysis of the IPA software (Qiagen, https://www.qiagenbio-informatics.com/products/ingenuity-pathway-analysis).

To generate the CRISPR-KO, single guide RNAs flanking the −50 kb Klf2 enhancer were designed by University of Michigan Transgenic Core (5′-CTACTACTTGGCAGGTTGGAGGG-3′ and 5′-GTCAAAGGGACCTGGTAGTTTGG-3′). Guide RNAs were tested for inducing chromosome breaks before microinjection. We screened 114 potential founders with PCR primers spanning the deletion (5′-ATGTGTGTGCATCTGGGGAGCAGAG-3′ and 5′-CCAGAGTGACTTTTCAGGCACAGGGG-3′) which generated a 450 bp product for the deleted allele. Primers within the deleted region were used to confirm a true indel (WT 5′-CTTATAACCTCCATTTCCTCCTCTGGG-′3 and WT 3′-CTTCGTGGTTTCCTGCTTGCTAAGATG-′3) that generates a 350 bp product for the WT allele. PCR products from 31 positive founders were cloned and sequence verified to characterize the deletion. A probe for Southern blot was designed using the following primers to clone out a 332 bp fragment from murine genomic DNA: 5′-CAAGGCCTTCCAGTACCAGG-3′ and 5′-TCTCAGTGGAGCTTGCTGTG-3′. The probe detects an RFLP in EcoRV-digested genomic DNA, 9.5 kb in WT allele and 5.6 kb in the CRISPR deleted allele [Klf2 ^{Δ-50kb(3.9kb)}]. Selected founders were outcrossed for two generations before being bred to homozygosity.

Mouse transgenesis at LBNL was performed in Mus musculus FVB strain mice. Animals of both sexes were used in these analyses and mouse embryos were excluded from further analysis if they did not encode the reporter transgene or if the developmental stage was not correct. For validation of in vivo enhancer activities, random Hsp68-lacZ transgenesis (for Klf2 elements) and enSERT was used for site-directed insertion of transgenic constructs at the H11 safe-harbor locus (Osterwalder et al., 2022; Kvon et al., 2020). EnSERT is based on pronuclear co-injection of Cas9, sgRNAs and a H11-homology arms-containing targeting vector encoding a candidate enhancer element upstream of a minimal promoter and a reporter protein (Kvon et al., 2020; Osterwalder et al., 2022). Related genomic enhancer coordinates are listed in Table 2. Predicted enhancer regions were PCR-amplified from mouse genomic DNA from WT FVB mice and cloned into a modified targeting vector encoding either an Hsp68-lacZ cassette (for random integration) or a human β-globin (HBB) minimal promoter upstream of a lacZ reporter (for enSERT). Embryos were excluded from further analysis if they did not contain a reporter transgene. CD-1 females served as pseudo-pregnant recipients for embryo transfer to produce transgenic embryos which were collected at E11.5 and stained with X-gal using standard techniques (Osterwalder et al., 2022). Embryos were harvested from timed matings at the timepoints indicated and pre-fixed in 2% paraformaldehyde-0.2% glutaraldehyde and stained as previously described (VanDusen et al., 2014a; Vincentz et al., 2019). After overnight staining at room temperature, embryos were post-fixed in 4% paraformaldehyde before imaging and sectioning. The number of tandems over total transgenics confirmed the negative activity of these elements as obtaining n=2 tandems (PCR-determined multicopy insertions at H11) for enSERT is sufficient to conclude whether an element is active or inactive (Kvon et al., 2020; Osterwalder et al., 2022).

If stained, β-galactosidase-stained embryos were post-fixed, washed in PBS, dehydrated, embedded, sectioned and NFR stained as previously described (George and Firulli, 2021; Vincentz et al., 2019). Images were acquired on the Keyence BZ-X800 florescence microscope system or the Leica DM5000 B compound microscope.

Conserved non-coding putative HAND2 binding regions were cloned out from genomic mouse DNA using the following primers: Igf2 5′-GAGAAGCTGGCAGATCAGGCTGTG-3′ and 5′-TGCTTCTGTTGAGAGGAGACAGTCTGG-3′; Igf2r 5′-TTGCCTGCATGTAAGTGTGCCTGG-3′ and 5′-TGTCTCTCAGGCTTCCTGTCTGGC-3′; Ptn 5′-ATTTCAGCTGGACTGCCATGGCAG-3′ and 5′-GGCTGGAAGAGGAGGCAAACAGAG-3′; Tmem108 5′-CATCATCACCATCACCATCGTCGTCG-3′ and 5′-GTATGCAGTGGACCTCTTTGACTTGTCAG-3′; Klf2 enhancer −16 kb element 5′-ATCTGTCCACCTCTACCTTCCA-3′ and 5′-AGTGGCTCTGACAACCTGAGAT-3′; Klf2 enhancer −50 kb element 5′-TGAACCTCCATTGATACACACC-3′ and 5′-GTCCCTAAGGATCATGTTGAGC-3′. Amplified sequences were Gibson (New England Biolabs) cloned into the pCR4-bG::lacZ-H11 enSERT vector and used to generate F0 enhancer transgenics. Briefly, the enSERT system uses CRISPR/Cas9-mediated site directed transgenesis at the murine H11 locus resulting in genomic integration of the human β-globin promoter with the enhancer element to be tested and the lacZ reporter cassette (Kvon et al., 2020). F0 embryos were harvested at E11.5 for β-galactosidase staining and analysis.

To generate the −50 kb Klf2 stable transgenic allele, primers corresponding to genomic region chr8:74791237-74793083 (mm9) (5′-AAGGGCCAGATGTGCTGAAA-3′ and 5′-GGCTGGTCTCGAACTCACAA-3′) were cloned into the HSP68-lacZ vector backbone as described previously (Vincentz et al., 2019) and used to create stable β-galactosidase-expressing mouse transgenic lines.

Section ISH was performed on 10 μm paraffin sections as described previously (George and Firulli, 2021). Whole mount ISH was performed using E10.5 embryos as described previously (George and Firulli, 2021). Antisense digoxygenin-labeled riboprobes were synthesized using T7, T3 or SP6 polymerases (Promega) and DIG-Labeling Mix (Roche) using the following plasmid templates: Hand2, Klf2.

Total RNA was isolated from E11.5 ventricles using the High Pure RNA Isolation Kit (Roche). RNA was used to synthesize cDNA using the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems). For qRT-PCR, cDNA was amplified using TaqMan Probe-Based Gene Expression Assays (Applied Biosystems) to quantify gene expression. qRT-PCR reactions were run on the QuantStudio 3 Real-Time PCR System (Thermo Fisher Scientific). Normalization to glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used to determine relative gene expression and statistical analysis was automatically applied by the instrument software. Significance of qRT-PCR results were determined by an unpaired two-tailed Student's t-test followed by post hoc Benjamini–Hochberg FDR correction as automatically calculated by the QuantStudio3 qRT-PCR thermal cycler software analysis package. Data are presented as relative quantitation values, where error bars depict the maximum and minimum values of each series of samples. A minimum n of 8 is used in all assays.

For ChIP assays, NIH3T3 cells were transfected with Lipofectamine 3000 with plus reagent (Invitrogen) according to the manufacturer's instructions with pCS2+Myc-Hand2, pCS2+Myc-E12 or pCS2 control constructs as indicated. After culturing for 48 h, SimpleChIP plus enzymatic chromatin IP kit (Cell Signaling Technology) was used as per the manufacturer's recommendations, and PCR was used to detect ChIP products run out on agarose gel.

We thank Danny Carney and Chloe Ferguson for technical assistance. The Klf2 expression vector used to generate the ISH probe was a kind gift from Jonathan A. Epstein. We thank Fabrice Darbellay for the generation of the β-globin-lacZ H11-targeting vector backbone and Nathan VanDusen for helpful comments.