Wt1 transcription factor impairs cardiomyocyte specification and drives a phenotypic switch from myocardium to epicardium

The heart is one of the first organs to acquire its function and it starts beating long before cardiac development is completed. In mammals, its function is essential to promote blood flow, and to sustain oxygenation and nutrition of the organism. Indeed, heart defects are among the major congenital anomalies responsible for neonatal mortality (van der Linde et al., 2011; Fitzgerald et al., 2015).

The zebrafish is a well-established vertebrate model organism in cardiovascular research, given its transparency during early developmental stages, its amenability to in vivo imaging and its rapid embryonic development (Eisen, 2020). Cardiac precursor cells derive from the anterior lateral plate mesoderm (Stainier et al., 1993). At 14 h post-fertilization (hpf), cardiac precursor cells start to express myosin light chain 7 (myl7) (Yelon et al., 1999) and sarcomere assembly begins soon after (Huang et al., 2009; Yang et al., 2014). As the assembly of sarcomeres continues, cardiac precursor cells migrate and fuse into a cone that later forms the heart tube, which is contractile at 24 hpf, and comprises a monolayer of cardiomyocytes lined in the interior with an endocardial sheet facing the lumen. Next, the heart tube starts to loop, leading to the formation of the two chambers: the atrium and the ventricle (Stainier et al., 1993). Concomitantly, more progenitors enter the heart tube through the arterial and venous poles (Knight and Yelon, 2016). Around 55 hpf, the outermost cell layer of the heart, the epicardium, starts to form. Epicardial cells arise from the proepicardium, a cell cluster derived from the dorsal pericardium that lies close to the venous pole of the heart. Cells from this cluster are later released into the pericardial cavity and attach to the myocardial surface, forming the epicardium (Peralta et al., 2013; Serluca, 2008).

Wilms tumor 1 (Wt1) is one of the main epicardial and proepicardial marker genes and plays a central role in epicardium morphogenesis (Moore et al., 1999; Serluca, 2008). Wt1 contains 4 DNA binding zinc-finger domains in the C-terminus and has been shown to act as a transcription factor (Hastie, 2017). Wt1 is transiently expressed in the epicardium during embryonic development and, in the adult heart, is reactivated after cardiac injury (van Wijk et al., 2012; Zhou et al., 2011).

The zebrafish has two Wt1 orthologues: wt1a and wt1b (Bollig et al., 2006). These genes are also expressed in the proepicardium and epicardium (Peralta et al., 2013; Serluca and Fishman, 2001) in partially overlapping expression domains. Using transgenic reporter and enhancer trap lines (Bollig et al., 2009; Perner et al., 2007), we have previously shown that wt1a and wt1b are initially expressed in a few proepicardial cells and later in epicardial cells (Andrés-Delgado et al., 2020; Peralta et al., 2014). Although wt1a and wt1b mRNA expression was not detected in the myocardium, wt1b:eGFP signal was transiently detected in cardiomyocytes of the atrium close to the inflow tract of the heart. Furthermore, some wt1a-associated regulatory regions were found to drive eGFP expression in cardiomyocytes (Peralta et al., 2014). Given that wt1a and wt1b regulatory elements drive gene expression in the myocardium, but endogenous mRNA expression is observed only in the proepicardium and epicardium, we hypothesized that wt1a and wt1b expression in the myocardium needs to be actively repressed to enable progression of normal heart development.

To explore whether there is a requirement for Wt1 ortholog downregulation in the myocardium for proper vertebrate embryonic development, we generated transgenic zebrafish models for tissue-specific overexpression of wt1a or wt1b in cardiomyocytes. We found that sustained wt1a or wt1b overexpression in the myocardium induced delamination and a phenotypic change from cardiomyocytes to epicardial-like cells. Moreover, we observed impaired cardiac morphogenesis, altered sarcomere assembly and delayed myocardial differentiation.

ATAC-seq data analysis of cardiomyocytes overexpressing wt1b revealed that this regulator acts as a break for cardiomyocyte differentiation by reducing chromatin accessibility of genomic loci associated with key processes such as sarcomere assembly, establishment of apicobasal polarity and adherens junctions formation. Altogether, our results demonstrate that transcriptional downregulation of wt1a and wt1b expression in cardiomyocytes is a prerequisite for cardiomyocyte specification, ensuring correct development of the heart and preventing a phenotypic switch of cardiomyocytes into epicardial cells.

During heart tube growth, cells from the pericardial mesoderm enter the heart tube at the venous pole (Knight and Yelon, 2016). These cells can be labeled with the line epi:eGFP (Peralta et al., 2013), an enhancer trap line of wt1a (Fig. 1A). We found that, during this process, cardiomyocyte precursors downregulate eGFP expression, concomitant with the activation of myl7:mRFP (Fig. 1B and Movie 1). We measured the eGFP/mRFP signal intensity ratio in cells of the heart tube, from the sinus venosus (SV) towards the growing heart tube. We found that the further away the cells were from the SV, the lower was the eGFP/mRFP ratio (Fig. 1C, n=3).

To further confirm our observations, we performed SMARTer RNA-seq of cells collected from three distinct regions: dorsal pericardium (dp), proepicardium (PE) and heart tube (ht) at 60 hpf (Fig. 1D-F). We detected a gradual decrease in wt1a and wt1b normalized counts among these three tissues, with the highest counts in PE cells and lowest in cells from the heart tube. The opposite trend was observed for myl7 expression, being highest in heart tube and lowest in the proepicardium samples (Fig. 1F).

Next, we performed lineage-tracing experiments for wt1a to see whether a pool of cardiomyocytes was derived from wt1a-expressing progenitors. We crossed the transgenic fish lines Tg(wt1a:CreERT2) with Tg(-3.5ubi:loxP-EGFP-loxP-mCherry) and induced recombination between 1 and 3.5 days post-fertilization (dpf). At 5 dpf, we analyzed which heart cells expressed GFP and are thus derived from wt1a-expressing progenitors. As expected, we found that most of the wt1a-derived cells were epicardial cells (Fig. 1G), but we also found a smaller population of cardiomyocytes that were derived from wt1a progenitors (Fig. 1G-H‴). In line with these results, previous work (Díaz del Moral et al., 2021) also identified a Wt1-derived cardiomyocyte pool in mouse hearts. Furthermore, based on single cell transcriptome analysis of embryonic mouse hearts, a common early progenitor cardiac pool for epicardial cells and cardiomyocytes was recently described (Tyser et al., 2021). We decided to study the temporal expression patterns of the epicardial marker genes Wt1 and Tcf21 in these mouse progenitor populations (Fig. S1A-A″). We found that Wt1 but not Tcf21 is expressed in the epicardial and myocardial common progenitor pool (Me5) at earlier stages. Like our lineage-tracing experiments, we observed that at later stages expression of Wt1 can also be found in a few cells in differentiated cardiomyocytes (Me3), suggesting that they are derived from Wt1 expressing Me5 progenitor cells. Moreover, during developmental stages, expression of Wt1 in Me5 decreases, indicating differentiation of the progenitors. These data suggest the existence of a common precursor population both in zebrafish and in mouse.

The observed downregulation of Wt1 expression during cardiomyocyte differentiation might indicate repression of Wt1 expression in mesodermal precursor cells during myocardial fate acquisition, while high levels are maintained in precursors destined to become epicardium. In line with this possibility there are regulatory elements of the wt1a locus that drive expression in the myocardium (Peralta et al., 2014). We wanted to further assess whether the observed downregulation of Wt1 during cardiomyocyte differentiation could be associated with directed repression of Wt1 gene expression. For this, we inspected previously published data on activating and repressing histone marks at the Wt1 genomic locus during four stages of cardiac differentiation from mouse embryonic stem cells (mESCs) (Wamstad et al., 2012) (Fig. S1B). Low levels of Wt1 transcripts were visible in mESCs but are absent throughout differentiation stages. Similarly, histone 3 K27 acetylation (H3K27ac) enrichment – which correlates with active promoter and enhancer activity – was present in regions proximal to the Wt1 transcriptional start site (TSS) in mESCs and mesodermal progenitor cells (Fig. S1B). Interestingly, some of these regions also colocalized with enhancer elements known to drive epicardial-specific reporter gene expression (Vieira et al., 2017). Conversely, histone H3 K27 trimethylation (H3K27me3) signatures, which associate with repressed regions, were nearly absent in mESCs, while in cardiac precursors cells and differentiated cardiomyocytes they massively decorated the extended regions flanking the Wt1 TSS, including the epicardial enhancer elements. Together, these observations indicate that during cardiomyocyte differentiation, Wt1 expression and epicardial enhancers become actively repressed. In summary, during heart development, cardiac precursor cells downregulate Wt1 orthologs upon their entry into the heart tube, which might be a prerequisite for their differentiation into cardiomyocytes (Fig. 1I).

We next aimed at exploring the biological relevance of the observed downregulation of Wt1 in cardiomyocytes. For this, we generated the line Tg(b-actin2:loxP-DsRED-loxP-eGFP-T2A-wt1a) to conditionally induce the expression of wt1a. Crossing this line into Tg(myl7:CreERT2) (Kikuchi et al., 2010) allowed the temporally induced overexpression of wt1a in cardiomyocytes. Hereafter, the double transgenic line is called myl7:CreERT2;eGFP-T2A-wt1a. We administered 4-hydroxytamoxifen (4-OHT) from 1 to 4 days post-fertilization (dpf) to induce recombination of loxP sites and activation of wt1a and eGFP expression during embryogenesis in cardiomyocytes (Fig. S2A). We confirmed wt1a and eGFP overexpression in the heart by RT-qPCR (Fig. S2B-E). Comparison of eGFP and wt1a expression between myl7:CreERT2;eGFP-T2A-wt1a with and without 4-OHT administration revealed a fourfold increase in eGFP and wt1a expression when 4-OHT was administered (Fig. S2B-E). Moreover, we generated a line to overexpress wt1b. We decided to use the Gal4/UAS system in this case, to allow a more homogeneous expression in the myocardium. By crossing the newly generated Tg(eGFP:UAS:wt1b) into Tg(myl7:Gal4) (Mickoleit et al., 2014), we specifically overexpressed wt1b and eGFP, under a bidirectional UAS promoter, in cardiomyocytes. The double transgenic line will be hereafter called myl7:Gal4;eGFP:UAS:wt1b (Fig. S2F). RT-qPCR data showed that wt1b and eGFP expression in the heart of the double transgenic line myl7:Gal4;eGFP:UAS:wt1b was upregulated fourfold compared with that of the single transgenic eGFP:UAS:wt1b (Fig. S2G,H). As a control, we used the double transgenic line Tg(eGFP:UAS:RFP);(myl7:Gal4) (Sanz-Morejon et al., 2019), hereafter named myl7:Gal4;eGFP:UAS:RFP. RT-qPCR analysis also indicated that expression of wt1a and wt1b could be monitored via GFP imaging (Fig. S2F-H).

Using these new lines, we analyzed the effect of sustained wt1a and wt1b overexpression in cardiomyocytes during heart development (Fig. 2A). In wt1a-overexpressing hearts but not in controls, we observed eGFP-positive cardiomyocytes located apically, protruding towards the pericardial cavity at 5 dpf (Fig. 2B-C‴). Moreover, these delaminating cardiomyocytes showed reduced expression of myosin heavy chain (MHC), suggesting the loss, to some extent, of the myocardial phenotype (Fig. 2C-C‴). We quantified how many of the delaminating cells were GFP⁺ or GFP⁻, and found that only GFP⁺ cells were delaminating, indicating that this delamination process is due to a cell-autonomous effect of wt1a in cardiomyocytes (Table S1). We detected a similar occurrence in wt1b-overexpressing hearts. Here, we also observed GFP-positive and MHC-negative cells that delaminated and adhered to the apical surface of the myocardium, starting at 3 dpf (Fig. 2D-E‴).

To better understand the origin of these apically positioned eGFP-positive cells in the wt1b overexpression hearts, we performed in vivo imaging in myl7:Gal4;eGFP:UAS:RFP and myl7:Gal4;eGFP:UAS:wt1b between 2 and 3 dpf (Fig. 2F,G and Movie 2). In myl7:Gal4:eGFP:UAS:wt1b hearts, some eGFP-positive cells started to round up and initiated delamination from the myocardium. Cells gradually changed from a flat to a rounded shape and ultimately remained adherent to the outer myocardial layer (Fig. 2G and Movie 2; n=4). This event of cell delamination was not observed in myl7:Gal4:eGFP:UAS:RFP control embryos (Fig. 2F and Movie 2; n=2). Apical extrusion of cardiomyocytes can be a consequence of myocardial malformation during which extruded cells are eliminated (Gentile et al., 2021; Rasouli et al., 2018). However, here we found that the delaminated cells remained attached to the myocardial surface. Between 5 and 6 dpf, these delaminated cells lost their rounded shape and flattened, acquiring an epicardial-like morphology (Fig. 2H,I and Movie 3, n=4).

To confirm that this type of cellular delamination with loss of MHC expression was specific to the overexpression of both wt1a and wt1b, we generated the Tg(eGFP:UAS:tcf21) to overexpress another well-known epicardial marker, tcf21 (Kikuchi et al., 2011), in cardiomyocytes (Fig. S3A,B,B′). We then crossed the Tg(eGFP:UAS:tcf21) into the Tg(myl7:Gal4) (Mickoleit et al., 2014) to overexpress tcf21 in cardiomyocytes. Contrary to what we observed when overexpressing wt1b in cardiomyocytes (Fig. S3C-E″), in the large majority of myl7:Gal4; eGFP:UAS:tcf21 embryos (59/62) we did not observe apical delamination. In the few cases where delamination occurred (3/62), the protruding cells still expressed MHC (Fig. S3D-D″). Moreover, to assess whether the detected changes were specific to the overexpression of wt1a and wt1b in cardiomyocytes, we overexpressed wt1b in other cardiac cells. We crossed Tg(eGFP:UAS:wt1b) with TgBAC(nfatc1:GAL4ff) (Pestel et al., 2016) to overexpress wt1b in the atrioventricular valves (Fig. S3F-F″), and with Tg(fli1a:Gal4) (Herwig et al., 2011) to overexpress wt1b in the endocardial cells (Fig. S3G-H″). We could not detect any apical delamination, looping defects or reduced MHC expression in these hearts at 3 dpf and 5 dpf.

Owing to the position and change of morphology of the delaminated cells, we hypothesized that these cells had undergone a change of fate. For a better characterization of a possible switch to an epicardial fate, we performed immunofluorescence labeling with the epicardial markers aldehyde dehydrogenase 2 (Aldh1a2) (Niederreither et al., 2002; Sugimoto et al., 2017) and caveolin 1 (Cav1) (Grivas et al., 2020) (Fig. 3A and Fig. S4A-I‴). We detected GFP/Aldh1a2 double-positive cells in wt1b overexpression hearts (n=3) but not in controls (n=4) (Fig. 3B-C″). Similarly, whereas in control hearts (n=6) we could not observe GFP/Cav1 double-positive cells (Fig. 3D-D″), in wt1b-overexpressing hearts (n=4), we identified GFP-positive cells that also expressed caveolin 1 (Fig. 3E-E‴). We also detected eGFP⁺ cells within the epicardium of wt1a-overexpressing hearts, but not in controls (Fig. 3F-I‴, and Fig. S4B-E‴). These eGFP⁺ cells did not express MHC and were Aldh1a2 positive (Fig. 3F-G‴ and Fig. S4D-E‴) as well as caveolin 1 positive (Fig. 3H-I‴ and Fig. S4F-I‴). We also tested for the colocalization of wt1a expression in cardiomyocytes with a third epicardial marker, transglutaminase b (tgm2b) (Weinberger et al., 2020). We performed in situ hybridization against tgm2b mRNA followed by immunohistochemistry against eGFP (Fig. S4J-M′). In non-recombined myl7:CreERT2;eGFP-T2A-wt1a hearts, tgm2b expression was visible in only a few epicardial cells and we could not observe any colocalization with eGFP-expressing cells (Fig. S4J-L′). However, in embryonically recombined myl7:CreERT2;eGFP-T2A-wt1a hearts, we could observe cells co-expressing tgm2b and eGFP located within the epicardium (Fig. S4K-M′).

These results suggest that, upon sustained ectopic overexpression of wt1a and wt1b, cardiomyocytes can delaminate apically from the myocardial layer and adopt features of epicardial cells that contribute to the formation of the epicardium.

We decided to get a better understanding on the cellular mechanisms underlying cardiomyocyte apical delamination upon wt1 overexpression (Fig. 4A). Previous reports showed that correct development and morphogenesis of the heart requires cell-cell adhesion and polarization of the cardiomyocytes (Phillips et al., 2007). The proper localization of tight junctions and adherens junctions has conventionally been used to assess the cell polarity (Zihni et al., 2016). We first performed immunostaining against ZO-1, a component of the tight junctions (Stevenson et al., 1986) (Fig. 4B-E′). Whereas the myl7:Gal4:eGFP:UAS:RFP control hearts (n=6) showed discrete apical localization of ZO-1 (Fig. 4B-C′), in myl7:Gal4:eGFP:UAS:wt1b hearts (n=8), ZO-1 levels were reduced, the signal was diffuse and not clearly localized to apical junctions between cardiomyocytes (Fig. 4D-E′). This suggests defects in the formation and localization of tight junctions upon wt1b overexpression. To evaluate the formation of adherens junctions, we crossed the Tg(myl7:cdh2-tdTomato) line (Fukuda et al., 2017) with myl7:Gal4:eGFP:UAS:wt1b. This allowed us to specifically visualize subcellular localization of cdh2-tdTomato in wt1b-overexpressing cardiomyocytes and control siblings (Fig. 4A,F-M′). At 5 dpf, in control embryos (n=6), tdTomato signal was clearly localized to cell-cell junctions (Fig. 4F,G) and detected apically in cardiomyocytes (Fig. 4H-I′). In contrast, myl7:cdh2-tdTomato;myl7:Gal4:eGFP:UAS:wt1b hearts (n=8) showed a diffuse and patchy staining for cdh2-tdTomato that was not restricted to the apical side of the cardiomyocytes (Fig. 4J-J″,L-M′). Moreover, we observed loss of cdh2-tdTomato signal in the delaminating cells, further indicating a loss of polarity in these extruding cells (Fig. 4K-K″). To confirm the impairment in the formation of adherens junctions, we carried out immunostaining against β-catenin (Fig. 4N-Q′), a core component of adherens junctions (Sheikh et al., 2009). Similar to what we had observed for cdh2-tdTomato, β-catenin staining was located at the apical side of cardiomyocytes in myl7:Gal4:eGFP:UAS:RFP control hearts (n=5) (Fig. 4N-O′). However, in wt1b-overexpressing hearts (n=5), β-catenin staining was no longer detected (Fig. 4P-Q′). Taken together, these data show that sustained expression of wt1b in cardiomyocytes leads to the mislocalization of tight junctions and adherens junctions, indicating an impairment of the apical domain in cardiomyocytes.

To understand the basal domain landscape of cardiomyocytes, we performed immunostaining against laminin (Fig. 4R-U′), a component of the basement membrane. Laminins have been associated with myocardial differentiation and with regulating the sarcolemmal properties (Derrick et al., 2021; Oliviéro et al., 2000; Wang et al., 2019, 2006; Yarnitzky and Volk, 1995). At 5 dpf, in the hearts of control fish (n=5), we observed clear anti-laminin staining at the basal and lateral domains of cardiomyocytes (Fig. 4R-S′), which correlates with previous observations (Derrick et al., 2021; Oliviéro et al., 2000). Laminin expression levels were severely reduced in wt1b-overexpression hearts (n=5), with no laminin observed in the lateral domains of the cardiomyocytes (Fig. 4T-U′). Thus, the observed reduced levels of laminin and its impaired deposition upon wt1b overexpression point towards an improper basal domain of cardiomyocytes. Taken together, our observations indicate that cardiomyocyte apicobasal polarization may be disrupted upon wt1b overexpression.

The disruptions in cell junctions and cell extrusion that we observed in wt1b-overexpressing cardiomyocytes led us to question the maturation and general architecture of these cells. Using whole-mount immunofluorescence, we observed reduced MHC staining in wt1b-overexpressing hearts at 1 dpf (n=6), when compared with controls (n=4) (Fig. 5A-C′). The reduction of MHC staining was specific to the heart, as it was not observed in the skeletal muscle of the myotome (n=6) (Fig. 5D). Although at 6 dpf we observed an increase in the levels of MHC signal in wt1b-overexpressing cardiomyocytes, the levels never reached those observed in the control group (Fig. 5E-G). We also analyzed myl7 mRNA expression levels using whole-mount in situ hybridization. Consistent with the results obtained using MHC immunostaining, at 3 dpf, myl7 expression was reduced in myl7:Gal4;eGFP:UAS:wt1b (23/25) compared with their single transgenic eGFP:UAS:wt1b control siblings (Fig. 5H,I). We reasoned that the reduced levels in MHC and myl7 staining could be indicative of an impaired maturation of cardiomyocytes. To test this hypothesis, we performed immunofluorescence staining against Alcam, a marker for undifferentiated cardiomyocytes (Hirata et al., 2006; Valente et al., 2019). At 6 dpf, we observed higher Alcam staining levels in hearts overexpressing wt1b (n=7) when compared with control hearts (n=6) (Fig. 5J-L).

We next analyzed whether sarcomere assembly was impaired in myl7:Gal4;eGFP:UAS:wt1b animals. We performed immunofluorescence staining against actinin, a protein known to be produced in the z line of the sarcomeres (Costa et al., 2002). Qualitative assessment of α-actinin revealed that not only were the levels lower but also the z-lines were thicker and shorter (Fig. 5M-P′) upon myocardial wt1b overexpression. Z-line disruption was particularly evident in delaminating cardiomyocytes (Fig. 5O,O′). We next sought to analyze sarcomere structure more in detail using serial block face scanning electron microscopy (SBFSEM) (Fig. 5Q-R′ and Movies 4-7). Z-bands were present at the sarcomere boundaries in both groups. While sarcomeres could be easily followed from z-band to z-band in the control heart (Fig. 5Q,Q′ and Movie 5), this was not possible in wt1b-overexpressing hearts (Fig. 5R,R′ and Movie 7). A further ultrastructural defect we observed in the wt1b overexpression heart was the presence of large intercellular spaces of extracellular matrix, between cardiomyocytes, the epicardium and endocardium. Moreover, while the control heart revealed a clearly visible basement membrane between the epicardium and the myocardium, as well as between the endocardium and myocardium (black line), this structure was not always visible in the wt1b-overexpressing heart (Fig. 5Q,R and Movies 4 and 6). This observation correlates with the impairment in laminin staining reported in wt1b-overexpressing cardiomyocytes (Fig. 4R-U′).

As we saw that cardiomyocyte structure and maturation were disrupted, we decided to evaluate cardiac performance. We performed in vivo imaging and analyzed different parameters for heart function in myl7:Gal4;eGFP:UAS:wt1b and myl7:Gal4; eGFP:UAS:RFP larvae (Fig. S5; n=14). myl7:Gal4;eGFP:UAS:wt1b ventricles presented a reduced stroke volume at 2 dpf (0.11±0.04 nl versus 0.04±0.03 nl) and this impairment did not recover at 5 dpf (0.39±0.17 nl versus 0.22±0.08 nl) (Fig. S5B). We next analyzed the heart rate. Although at 2 dpf we could not detect changes in heart rate [114±8 beats per min (bpm) versus 119±8 bpm] we observed a significant decrease in heart frequency at 5 dpf upon wt1b overexpression (166±13 bpm versus 141±9 bpm) (Fig. S5C). At 2 dpf, the ejection fraction did not significantly change between either group (43±14% versus 51±8%). However, at 5 dpf there was a clear reduction in the ejection fraction of the atrium in wt1b-overexpressing hearts compared with controls (55±8% versus 41±12%). In contrast, the ventricular ejection fraction was initially significantly reduced at 2 dpf in wt1b-overexpressing embryos (48±13% versus 35±13%) but recovered at 5 dpf (50±8% versus 49±8%) (Fig. S5D). The loss of cardiac function suggests that wt1b overexpression led to heart failure because, on average, only 28% of myl7:Gal4;eGFP-UAS-wt1b larvae survived past 6 dpf (of a total of 600 embryos from six independent experiments).

Altogether, our findings indicate that sustained expression of wt1b in cardiomyocytes impaired cardiomyocyte maturation, including cardiomyocyte sarcomere assembly and the extracellular matrix. This ultimately disrupted heart function and reduced survival.

Seeing that wt1b overexpression in cardiomyocytes induced several cardiac malformations and caused a phenotypic change in some cells, we decided to explore how the sustained expression of this transcription factor affected chromatin accessibility. To achieve this, we performed an assay for transposase-accessible chromatin sequencing (ATAC-seq) (Buenrostro et al., 2013) in 5 dpf, FACS-sorted GFP⁺ cells from either the myl7:Gal4;eGFP:UAS:RFP control or myl7:Gal4;eGFP:UAS:wt1b larvae (Fig. 6A). We identified 1452 differential peaks in wt1b-overexpressing cardiomyocytes, of which almost all, except for 14 peaks, showed reduced chromatin accessibility (Fig. 6B and Table S2). Most of the differential accessible regions were located close to promoter regions (38.87%), in introns (30.37%) or in distal intergenic regions (26.14%) (Fig. 6C). We performed gene ontology (GO) analysis for the genes lying near the differentially accessible regions. From the top 25 biological pathways that had reduced accessibility of regions close to genes associated with these pathways (Fig. S6), five of them account for muscle development (Fig. 6D and Table S3). Of the top 25 cellular component pathways (Fig. S6) we found some to be involved in ‘actin cytoskeleton’, ‘basolateral plasma membrane’, ‘apical part of the cells’, ‘contractile fiber’ or ‘myofibril’ (Fig. 6E and Table S3). From the top 25 molecular function pathways (Fig. S6) five of them are directly implicated in transcription regulation and another four in cytoskeleton formation and cell adhesion, such as ‘actin binding’, ‘actin-filament binding’, ‘cell adhesion molecule binding’ and ‘beta-catenin binding’ (Fig. 6F and Table S3). All of these pathways, which are underrepresented in the myl7:Gal4;eGFP:UAS:wt1b samples, strongly correlate with the defects observed in hearts overexpressing wt1b. To identify potential transcription factors that might be binding to the differentially accessible regions, we performed MEME-Centrimo motif analysis, and found Wt1 to be one of the top 5 motifs represented (E-value=7.8e-5) (Fig. 6G). This motif could be identified in 672 (46.25%) of the differentially accessible regions (three open and 669 closed) (Fig. 6G′). To further investigate which of the open regulatory regions and their associated genes were potential direct targets of Wt1, we compared our ATAC-seq data with Wt1 target genes identified in the CHIP-atlas database (Oki et al., 2018). 41% of the regions associated with differential accessibility identified in our ATAC-seq (426, of which only six represent regions with open chromatin) were shared with the CHIP-atlas database for Wt1. GO analysis of the associated common genes identified pathways similar to those observed previously, suggesting a direct regulation of these pathways by Wt1b (Fig. 6H).

Having seen that overexpression of wt1b in cardiomyocytes affected heart development and that these changes correlated with the observed molecular signature, we looked more closely at how the genetic landscape of some of the genes, with associated differentially accessible regions, was affected. We had previously seen that apical cell-cell junctions were disrupted in the hearts of the myl7:Gal4;eGFP:UAS:wt1b embryos, including expression and localization of cdh2, ZO-1 and β-catenin. In agreement, we observed that putative regulatory regions near cdh2 and ctnna (another core component of adherens junctions) revealed lower accessibility when wt1b was overexpressed in cardiomyocytes (Fig. 6I,I′). We also observed lower accessibility in core apicobasal polarity pathway genes (Assémat et al., 2008) such as pard6b and pard3bb from the apical polarity pathways, and scrib, dlg1a and dlg1b from the basolateral pathways (Fig. 6J,K and Table S4), supporting a perturbed apicobasal polarity in the wt1b overexpression lines. Moreover, we detected that several genes associated with sarcomere assembly, such as e2f3, rbfox2, rybpa and rybpb (Gallagher et al., 2011; Henry et al., 2020; King et al., 2008), presented lower chromatin accessibility in wt1b-overexpressing cells (Fig. 6L-L′ and Table S4), which could explain the disrupted sarcomeres observed (Fig. 5M-R).

In conclusion, ATAC-seq data analysis revealed that wt1b overexpression in the heart decreased overall chromatin accessibility associated with key genes involved in cardiomyocyte maturation and structural differentiation, with Wt1b likely to directly repress gene expression programs controlling muscle development, cell polarity and actin binding.

During myocardial development, cells from the precardiac mesoderm enter the heart tube and contribute to its growth. We found that myocardial expansion at the venous pole occurs through a population of wt1a and wt1b-positive cells that downregulate wt1a or wt1b expression during differentiation. We found that sustained wt1b activity in cardiomyocytes reduced chromatin accessibility in regulatory regions associated with cardiomyocyte-specific genes and that this can induce their phenotypic switch from myocardial to epicardial-like cells. In a different context, Wt1 expression has been shown to prevent the activation of a muscle differentiation program (Miyagawa et al., 1998). Wt1 also impairs the differentiation of embryonic stem cells towards a myocardial lineage (Wagner et al., 2021). These observations suggest that Wt1 orthologs might be generally repressing myocyte specification. However, Wt1 loss of function in cardiomyocytes also impairs heart development (Díaz del Moral et al., 2021) underlining the need for a delicate balance of Wt1 levels in cardiomyocytes during heart formation.

In this work we report the effect of the overexpression of the –KTS isoform of Wt1, that can act as a transcriptional activator or repressor (Toska and Roberts, 2014). Our ATAC-seq results suggest that in cardiac progenitors, Wt1 might act as a repressor of the myocardial gene program. It will be interesting to further analyze if the +KTS isoform, associated with RNA-binding properties (Toska and Roberts, 2014) would have similar effects on myocardial specification.

While a subset of pericardial cells enters the heart tube, downregulate wt1a and wt1b, and contribute to myocardial expansion, other cells maintain the expression and form the proepicardium. During proepicardium formation, wt1a and wt1b-positive cells extrude apically from the dorsal pericardial mesothelium giving rise to proepicardial cell clusters that subsequently are transferred to the myocardium forming the epicardium (Andrés-Delgado et al., 2019). Here we find that myocardial wt1a or wt1b-positive cells undergo a similar process and delaminate apically from the myocardium. It will be important to further decipher possible parallelisms between these two processes and elucidate the direct role of wt1a and wt1b during these cellular rearrangements. Wt1 participates in the mesothelial-to-mesenchymal transition giving rise to epicardial-derived cells (EPDCs) (Martínez-Estrada et al., 2010; von Gise et al., 2011). Moreover, Wt1 has been suggested to control the retinoic acid (RA) signaling pathway during EPDC formation (Guadix et al., 2011; von Gise et al., 2011). The fact that cardiomyocytes overexpressing wt1a or wt1b are relocating to the epicardial layer might indicate that these cells undergo an EMT-like process that might be mediated by RA. However, the absence of aldh1a2 expression in the myocardium, prior to delamination suggests that aldh1a2 expression might be a consequence rather than a cause of apical delamination of wt1a or wt1b expressing cells. Of note, not all eGFP-positive cardiomyocytes undergo delamination. It might thus be possible that not all cardiomyocytes have the capacity to respond to the same extent to wt1a or wt1b overexpression. Indeed, in the mouse a small subset of cardiomyocytes has been shown to express Wt1 and as such, not all cardiomyocytes might be equally sensitive to a change in Wt1 dosage (Cano et al., 2016; Rudat and Kispert, 2012). Furthermore, this effect is very specific to wt1a and/or wt1b function in cardiomyocytes and has not been observed neither by expressing wt1b in other cardiac cell types nor by overexpressing tcf21 in cardiomyocytes. Cardiomyocyte extrusion has been observed in klf2 and snai1 mutant zebrafish (Gentile et al., 2021; Rasouli et al., 2018). While in these mutants extruded cardiomyocytes are eliminated, here we report that the extruded wt1a or wt1b-positive cells remain on the myocardial surface contributing to the epicardial layer.

Wt1 overexpression in cardiomyocytes had been suggested to trigger a change in cell fate in arrhythmogenic right ventricular cardiomyopathy (ARVC) (Dorn et al., 2018). In this disease, a subset of cardiomyocytes has been proposed to start to express Wt1 and convert into adipocytes. Interestingly, epicardial fat represents an epicardial derivative (Zangi et al., 2017). Together with our results, this indicates that expression of Wt1 in cardiomyocytes contributes to a phenotypic change, transforming them into epicardial cells or EPDC-like cells. Genetic profiling of these delaminated cells might provide a better insight into the changes that these cells undergo, to allow for their delamination and phenotypic switch. Previous findings also pointed to the possibility of cell plasticity between epicardium and myocardium whereby epicardial cells could transdifferentiate into cardiomyocytes during development and repair (Smart et al., 2011; Zhou et al., 2008). Here we present evidence of the opposite phenotypic switch, triggered by one single transcription factor. Our work shows that during cardiac development, wt1a or wt1b expression is turned off in cardiomyocytes once they enter the heart tube to allow their correct differentiation. Dissecting the regulatory mechanisms controlling Wt1 ortholog transcription in cardiomyocyte precursors will further expand our knowledge on the tight spatio-temporal control of heart tube expansion and concomitant cardiomyocyte differentiation.

Experiments were conducted with zebrafish (Danio rerio) embryos and adults aged 3-18 months, raised at maximal 5 fish/l. Fish were maintained under the same environmental conditions: 27.5-28°C, with 14 h of light and 10 h of dark, 650-700 μs/cm, pH 7.5 and 10% of water exchange daily. Experiments were conducted after the approval of the ‘Amt für Landwirtschaft und Natur’ from the Canton of Bern, Switzerland, under the licenses BE95/15 and BE 64/18. All transgenic lines used in this manuscript can be found in Table S5.

To generate the transgenic line eGFP:UAS:wt1b and the eGFP:UAS:tcf21, the RFP fragment from the plasmid used to clone eGFP:5xUAS:RFP (Sanz-Morejon et al., 2019) was replaced by either the coding sequence of wt1b(-KTS) isoform or of tcf21, PCR amplified from 24 hpf and 5 dpf zebrafish embryo cDNA and assembled using Gibson cloning. The final entire construct is flanked with Tol2 sites to facilitate transgenesis. In this line, tissue-specific expression of Gal4 drives the bidirectional transactivation of the UAS leading to the expression of both eGFP and wt1b(-KTS) or tcf21-coding sequence. The full name of these lines is Tg(bGI-eGFP:5xUAS:wt1b(-KTS)-bGI;cryaa:eCFP)^brn4, Tg(bGI-eGFP:5xUAS:tcf21)-bGI; cryaa:eCFP)^brn5.

The construct bactin2:loxP-DsRed2-loxP-eGFP-T2A-wt1a was generated by Gateway cloning (MultiSite Gateway Three-Fragment Vector Construction Kit; Invitrogen). As a destination vector, pDestTol2pA2 was used. The floxed DsRed2 cassette was derived from vector pTol2-EF1alpha-DsRed(floxed)-eGFP (Hans et al., 2009) and the wt1a cDNA was amplified from vector pCS2P-wt1a (Bollig et al., 2006). The final construct is flanked with Tol2 sites to facilitate transgenesis. In the resulting zebrafish line, DsRed is expressed from the ubiquitous β-actin promoter. After Cre-mediated excision of the STOP cassette, both eGFP and wt1a are expressed in a tissue-specific manner. The full name of this line is Tg(bactin2:loxP-DsRed2-loxP-eGFP-T2A-wt1a)^li21.

4-hydroxytamoxifen (4-OHT; Sigma, H7904; Table S5) stock was prepared by dissolving the powder in ethanol to 10 mM concentration. To aid with the dissolution, the stock was heated for 10 min at 65°C and then stored at −20°C, protected from the light. 4-OHT was administered at the indicated times, at a final concentration of 10 µM. For embryos, treatments were performed continuously. Prior to administration, the 10 mM stock was warmed for 10 min at 65°C (Felker et al., 2016).

For in vivo imaging of the beating zebrafish heart, 2 dpf embryos were pipetted with melted 1% low melting agarose in E3 medium (about 45°C), containing 0.003% 1-phenyl-2-thiourea (PTU) (Sigma-Aldrich) to avoid pigmentation and tricaine at 0.08 mg/ml (pH 7) to anaesthetize the fish, into a U-shaped glass capillary (Leica). This U-shaped capillary was mounted in a 35 mm MatTek imaging dish. The dish was filled with E3 medium containing 0.003% PTU and tricaine at 0.08 mg/ml (pH 7).

Imaging was performed with the Leica TCS SP8 digital light sheet (DLS) microscope. We used a 25× detection objective with NA 0.95 water immersion and a 2.5× illumination objective with a light sheet thickness of 9.4 µm and length of 1197 µm. The total field of view is 295×295 µm, fitting the size of the embryonic zebrafish heart, allowing space for sample drift. The images were acquired in XYTZL acquisition (XY, single optical section; T, time series; Z, serial optical sections; L, looped acquisition) mode for later retrospective gating. The parameters as shown in Table S6 were applied.

The images were saved as single lif file and transferred to a workstation (HP-Z series, Dual Intel Xeon e5-2667 v4 3.2 GHz, 256 GB, NVIDIA GeForce GTX 1080 Ti). A quality check of the data was performed before the data were further processed. The survival of the larva until the end of the acquisition, the sample drift and the degree of bleaching were assessed in the Processor (https://github.com/MercaderLabAnatomy/PUB_Marques_et_al_2021/tree/main/ImageProcessing_BeatingHeartReconstruction). The data were only used if the larva survived the acquisition. The single lif file was converted to XYTC tif files, using the Converter_6D (https://github.com/MercaderLabAnatomy/PUB_Marques_et_al_2021/tree/main/ImageProcessing_BeatingHeartReconstruction). Each XYTC file was named in the following format ‘Image_R0000_Z0000’ to be recognized for further processing.

Retrospective gating was performed as previously described (Liebling et al., 2005, 2006; Ohn et al., 2009). The MATLAB (R2017a) tool BeatSync V2.1 was used for retrospective gating (access to the software can be requested from the research group of Michael Liebling, IDIAP Research Institute, Martigny, Switzerland). The settings for re-synchronization in the BeatSync software were ‘Normalized mutual information’, ‘Recursive Z-alignment’ and ‘Nearest-neighbor interpolation’. One entire heart cycle was re-synchronized in 3D. After re-synchronization, a 3D time lapse of a virtually still heart was created, using the Fiji (Schindelin et al., 2012) tool Make_timelapse (https://github.com/MercaderLabAnatomy/PUB_Marques_ et_al_2021/blob/main/ImageProcessing_BeatingHeartReconstruction/Make_timelapseMovie.py). The time lapse was represented as maximum intensity projection or individual optical slices.

Dorsal pericardium, proepicardium and heart tube were manually dissected, with tungsten needles, from 60 hpf epi:eGFP;myl7:mRFP zebrafish larvae. A minimum of 10 of each tissue/organ were collected for each sample in ice-cold PBS. Cells were centrifuged for 7 min at 250 g. The excess liquid was removed and the cells were stored at −80°C until further use. RNA was directly transformed and amplified into cDNA from the lysed tissue using the SMARTer Ultra Low Input RNA for Illumina Sequencing – HV kit. cDNA quality control was verified using the Agilent 2200 BioAnalyzer and the High Sensitivity DNA Chip from Agilent's High Sensitivity DNA Kit. Next, 50 ng of amplified cDNA were fragmented with the Covaris E220 (Covaris) and used for preparing sequencing libraries using the TruSeq RNA Sample Prep Kit v2 kit (Illumina), starting from the end repair step. Finally, the size of the libraries was checked using the Agilent 2200 Tape Station DNA 1000 and the concentration was determined using the Qubit fluorometer (Thermo Fisher Scientific). Libraries were sequenced on a HiSeq 2500 (Illumina) to generate 60 base single reads. FastQ files were obtained using CASAVA v1.8 software (Illumina). NGS experiments were performed in the Genomics Unit of the Centro Nacional de Investigaciones Cardiovasculares (CNIC).

All bioinformatics analyses were performed using bash scripts or R statistical software. A quality check of the samples was performed using FASTQC and reports summarized using MultiQC (Andrew, 2010; Ewels et al., 2016). Adapters from the fastq files were trimmed using fastp software (Chen et al., 2019). Reads were aligned to GRCz11 danRer11 v102 assembly from Ensembl using STAR (Yates et al., 2020). The reads were summarized using featureCounts (Liao et al., 2014). The counts data were imported to Deseq2 and genes who had expression across all samples (rowSums) greater than or equal to 10 were kept, ensuring the reliable expression estimates (Love et al., 2014). After evaluation of the PCA, one of the samples from the heart tube was determined as an outlier and removed from the downstream analysis. The differential expression analysis was performed using ‘ashr’ LFC Shrinkage (Stephens, 2017). A gene was considered as significant if the P adjusted value was less than 0.05. The plots were plotted using ggplot2 (Wickam et al., 2016).

Whole-mount immunofluorescence on embryos was carried out as previously described (Sanz-Morejon et al., 2019). In brief, embryos were fixed overnight at 4°C in 4% paraformaldehyde (PFA) (EMS, 15710), then washed with PBS-Tween20 (0.1%) and permeabilized for 30-60 min with PBS-TritonX100 (0.5%), depending on the stage and the antibody used. Permeabilization was followed by blocking for 2 h with histoblock (5% BSA, 5% goat serum and 20 mM MgCl₂ in PBS). Afterwards, embryos were incubated overnight at 4°C with the primary antibodies in 5% BSA. The next day, embryos were washed with PBS-Tween20 (0.1%) followed by an overnight incubation in the secondary antibodies at 4°C in 5% BSA. Finally, embryos were washed with PBS-Tween20 (0.1%) and a nuclear counterstaining with DAPI (Merck, 1246530100) 1:1000 was carried out.

Immunofluorescence on paraffin sections was performed as previously described (González-Rosa et al., 2011). Briefly, paraffin sections were dewaxed and rehydrated through a series of ethanol incubations, as previously described for histological staining. Afterwards, epitope was recovered by boiling the samples in 10 mM citrate buffer (pH 6) for 20 min. Next the same procedure was applied as described above for whole-mount immunofluorescence.

The following primary antibodies were used: anti-eGFP (Aves, eGFP-1010) at 1:300, anti-myosin heavy chain at 1:50 (DSHB Iowa Hybridoma Bank, MF20), anti-Aldh1a2 at 1:100 (Gene Tex), anti-alcama at 1:100 (DSHB Iowa Hybridoma Bank, Zn-8), anti-α-actinin at 1:200 (Sigma Aldrich), anti-caveolin 1 at 1:100 (BD Biosciences), anti-ZO1 at 1:200 (Invitrogen), anti-laminin at 1:20 (Sigma) and anti-β-Catenin at 1:100 (BD Biosciences). Secondary antibodies were Alexa Fluor 488, 568 and 647 (Life Technologies) at 1:250 and biotin anti-rabbit (Jackson Immuno Research, 111-066-003) followed by StreptavidinCy3 or Cy5 conjugate (Molecular Probes, SA1010 and SA1011) at 1:250. All catalog codes are in Table S5.

Hearts from Tg(eGFP:5xUAS:wt1bOE-KTS;myl7:Gal4) and Tg(eGFP:5xUAS:RFP;myl7:Gal4) were extracted at 40 dpf. Ventricle, atrium and bulbus arteriosus were manually dissected and stored separately in pools of five. For each sample, up to three biological replicates were collected. Tg(eGFP:UAS:tcf21; -1.5hsp70l:GAL4) embryos were heat shocked at 24 hpf and collected in pools of 50 embryos at 48 hpf for RNA extraction. Each pool of embryos represented a biological replicate. A total of four biological replicates was used for each condition. Total RNA was extracted by using TRI Reagent (Sigma-Aldrich; T9424) according to the manufacturer's recommendations. Afterwards, 200 ng of total RNA was reverse-transcribed into cDNA using High Capacity cDNA Archive Kit (Invitrogen Life Technologies; 4374966). Quantitative PCR (qPCR) was performed in a 7900HT fast real-time PCR system (Applied Biosystems). qPCR was carried out using Power Up SYBR Green Master Mix (Applied Biosystems, A25742). All primer details are in Table S5.

The PCR program was run as follows: initial denaturation step for 30 s at 95°C, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. To calculate the relative index of gene expression, we employed the 2^−ΔΔCt method, using e1f2a expression for normalization.

In situ hybridization on paraffin sections was carried out as follows: paraffin sections were dewaxed and rehydrated through a series of ethanol incubations. Sections were then refixed with 4% PFA at room temperature for 20 min. Afterwards, they were washed with PBS and the tissue was permeabilized by incubating the slides with 10 µg/ml of proteinase K, for 10 min. at 37°C. Afterwards, slides were washed with PBS and briefly refixed with 4% PFA. The tissue was then incubated for 10 min with triethanolamine 0.1 M (pH 8) and 0.25% acetic anhydride. After washing the slides with PBS and RNAse-free water, the slides were incubated for 3 h with pre-hybridization buffer [50% formamide, 5×SSC (pH 5.5), 0.1×Denhardt's, 0.1% Tween20, 0.1% CHAPs and 0.05 mg/ml tRNA] at 65°C. Afterwards, pre-hybridization buffer was replaced with hybridization buffer (pre-hybridization buffer with mRNA probe). Slides were left to incubate with hybridization buffer overnight at 65°C. The next day, slides were washed twice with post-hybridization buffer I [50% formamide, 5×SSC (pH 5.5) and 1% SDS] and twice with post-hybridization buffer II [50% formamide, 2×SSC (pH 5.5) and 0.2% SDS]. Each wash was carried out for 30 min at 65°C. Slides were then washed another three times with maleic acid buffer (MABT) and then incubated for 1 h blocking solution (2% fetal bovine serum, heat inactivated, and 1% blocking reagent, in MABT). Tissue was incubated overnight at 4°C with anti-DIG antibody in blocking solution at 1:2000. Finally, sections were thoroughly washed with MABT and incubated in alkaline phosphatase buffer [AP buffer, NaCl 0.1 M, MgCl₂ 0.05 M and 10% Tri-HCl (pH 9.5)]. Finally, colorimetric assay was performed using BM purple. After the desired staining was achieved, slides were washed with PBS and fixed with 4% PFA, before mounting them with 50% glycerol and imaging using a Zeiss Imager M2, with an Olympus UC50 camera.

After imaging, sections were washed and further permeabilized with PBS with 0.5% Triton X-100. They were then incubated for 2 h with 5% BSA at room temperature and incubated with primary antibody, chicken anti-GFP (1:300 in 5% BSA), overnight at 4°C. The next day, slides were washed in PBS-0.1% Tween20 and incubated for 1 h at room temperature with secondary antibody anti-chicken-HRP. Signal was obtained by incubating slides with DAB solution for 30 s at room temperature. The reaction was stopped with water. Slides were then mounted in 50% glycerol and imaged. Reagent details can be found in Table S5.

Whole-mount in situ hybridization was performed as described previously (Woltering et al., 2009), with some minor adaptations. Embryos were selected at 24 hpf and 3 dpf for eGFP expression. After fixation, the embryos were washed with PBS and gradually dehydrated through a methanol series. Embryos were stored in 100% methanol for a minimum of 2 h, at −20°C. Afterwards, the embryos were rehydrated and permeabilized with proteinase K (10 μg/ml in TBST) at 37°C. Incubation times were adjusted according to the stage of the embryos (24 hpf, 10 min and 72 hpf, 20 min). This was followed by a 20 min incubation in 0.1 M triethanolamine (pH 8) with 25 μl/ml acetic anhydride.

After 4 h of pre-hybridization at 68°C, myl7 riboprobe was diluted in pre-hybridization solution at a concentration of 300 ng/ml. The embryos were incubated with the riboprobe overnight at 68°C. The following day, the riboprobe was removed and the embryos were incubated twice for 30 min with post-hybridization solution at 68°C. Embryos were then incubated with blocking buffer, freshly prepared and afterwards incubated overnight with anti-DIG antibody (in blocking solution) at 1:4000 at 4°C.

The embryos were then washed extensively with maleic acid buffer [150 mM maleic acid (pH 7.5), 300 mM NaCl and 0.1% Tween 20]. Finally, the embryos were transferred to a 6-well plate and pre-incubated with AP-buffer [0.1 M Tris base (pH 9.5), 0.1 M NaCl, 1 mM MgCl2, 0.1% Tween 20] and then incubated with BM-purple at room temperature. As soon as color was visible in the heart of either group (overexpression or control), the staining was stopped in both groups by adding TBST and embryos were re-fixed in 4% PFA. Reagent details can be found in Table S5.

Using a microscope, we could obtain pictures of the hearts of the embryos. For image acquisition, embryos were mounted on 3% methylcellulose for ease of orientation. Embryos were positioned so that most of the heart could be observed in a single plane.

Images were acquired with a Nikon SMZ800N stereomicroscope. Illumination conditions and acquisition parameters were maintained for all embryos.

Zebrafish embryos at 5 dpf were euthanized with an 0.048% of tricaine and immediately fixed with 2.5% glutaraldehyde with 0.15 M cacodylate buffer and 2 mM CaCl₂ (pH 7.4). Embryos were then processed for serial block face scanning electron microscopy as previously described (Odriozola et al., 2017 preprint). Briefly, we proceed as follows: samples were rinsed three times in ice-cold 0.15 M Na-cacodylate for 5 min. They were then incubated in 0.15 M Na-cacodylate solution containing 2% OsO₄ and 1.5% potassium ferrocyanide for 45 min at room temperature and for 15 min in a water bath at 50°C. Samples were rinsed three times for 5 min in water. They were then incubated with 0.64 M pyrogallol for 15 min at room temperature, for 5 min in a water bath at 50°C and subsequently rinsed with water. The embryos were incubated in 2% OsO₄ for 22 min at room temperature and 8 min in a water bath at 50°C. Afterwards, they were rinsed again in water (three times for 5 min) and incubated overnight in a solution of 0.15 M gadolinium acetate (LFG Distribution) and 0.15 M samarium acetate (LFG Distribution; pH 7.0). The next day the embryos were rinsed three times for 5 min with water and incubated in 1% Walton's lead aspartate (Walton, 1979) at 60°C for 30 min and rinsed with water (3 times for 5 min).

After staining, the samples were dehydrated in a graded ethanol series (20%, 50%, 70%, 90%, 100% and 100%) at 4°C, with each step lasting 5 min. They were then infiltrated with Durcupan resin mixed with ethanol at ratios of 1:3 (v/v), 1:1 and 3:1, each step lasting 2 h. The embryos were left overnight to infiltrate with Durcupan. The next day, samples were transferred to fresh Durcupan and the resin was polymerized for 3 days at 60°C. Sample blocks were mounted on aluminum pins (Gatan) with a conductive epoxy glue (CW2400, Circuitworks). Care was taken to have osmicated material directly exposed at the block surface in contact with the glue to reduce specimen charging under the electron beam. Pyramids with a surface of ∼500×500 μm² were trimmed with a razor blade.

Three-dimensional (3D) ultrastructural images were produced by serial block face scanning electron microscopy (SBFSEM) on a Quanta FEG 250 SEM (FEI, Eindhoven, The Netherlands) equipped with a 3View2XP in situ ultramicrotome (Gatan). Images were acquired in low vacuum mode (40 mPa), except where indicated otherwise. Acceleration voltage was 5 kV and pixel dwell time was set to 2 µs. Image acquisition was carried out with a back scattered electron detector optimized for SBFSEM (Gatan). Image stacks were aligned, normalized and denoised by non-linear anisotropic diffusion in IMOD (Kremer et al., 1996). Each field of view consisted of 8192×8192 pixels with a dimension of 6 nm/pixel in x-y and 50-150 nm in the z direction. Final image montage was created in Fiji.

Heartbeat analysis was performed by assessing the following parameters: degree of rhythmic beating as root mean square of successive differences (RMSSD) (Collins et al., 2019); stroke volume (SV; difference between diastolic and systolic volume); ejection fraction (EF; difference between diastolic and systolic volume relative to the diastolic size); cardiac output (CO; SV multiplied by heart rate); and diastolic volume and heart rate as described previously (DeGroff, 2002; Yalcin et al., 2017).

We recorded 300 frames of the beating heart in the GFP channel in Tg(myl7:Gal4;eGFP:UAS:wt1b) and Tg(myl7:Gal4;eGFP:UAS:RFP) at 2 dpf and 5 dpf using the fluorescence stereo microscope Nikon SMZ25 (SHR Plan apo 1× objective, 10× zoom, 2880×2048 pixel, 0.44 µm/pixel, 17 frames/s).

We developed the FIJI plug-in ‘Heart beat analysis’ to sequentially process all images in a folder and guide the user through each manual step of the analysis. The manual steps are (1) find the two diastolic and systolic states of the heart, (2) adjust a line to D_L and D_S, and (3) draw one line at the border of the ventricle. The plug-in ‘Heart beat analysis’ opens subsets of the data (100 frames and only green channel per fish from the nd2 RGB file), applies a Gaussian blur filter (10 px), indicates which manual step to perform, calculates the HR by detecting maxima in a kymograph and subsequently saves all kymograph images as tiff files, results as csv files and all lines as zip files in ROI sets.

myl7:Gal;eGFP:UAS:RFP and myl7:Gal;eGFP:UAS:wt1b embryos at 5 dpf were used to obtain GFP⁺ heart cells. The heart region of these embryos was manually dissected and placed in Ringer's solution. Afterwards, the tissue was briefly centrifuged in a tabletop centrifuge and the Ringer's solution was replaced by a mix of 20 mg/ml collagenase in 0.05% trypsin. The samples were incubated at 32°C for 25 min. Every 5 min this mixture was gently mixed. The tissue was visually inspected for dissociation. After cell disaggregation, the reaction was stopped with Hanks' solution (1×HBS, 10 mM Hepes and 0.25% BSA). The homogenized samples were centrifuged at 250 g for 10 min and re-suspended in Hanks' solution. The cells were then passed through a 40 µm filter, centrifuged again for 10 min at 400 g and re-suspended in 50 µl of Hanks' solution for FAC sorting. Dead cells were marked with 7-aminoactinomycin D (Invitrogen) and discarded. Cells were FAC sorted into Hanks' solution, on a Moflo astrios EQ (Beckman Coulter), and analyzed for forward and side scatter, as well as for eGFP fluorescence. Between 1200 and 1500 cells per sample were sorted for ATAC-seq.

FAC-sorted GFP⁺ cells were gently centrifuged and Hanks' solution was replaced by lysis buffer [10 mM Tris-HCL (pH 7.4), 10 mM NaCl, 3 mM MgCl₂ and 0.1% IGEPAL CA-630]. Cells were immediately centrifuged at 500 g for 10 min at 4°C. The supernatant was discarded and replaced with the transposition reaction mix (Tn5 in TD buffer) for tagmentation, and incubated at 37°C for 30 min. Afterwards, 500 mM of EDTA was used for quenching. The solution was incubated for 30 min at 50°C. MgCl₂ was added to a final concentration of 45 mM. Samples were stored at 4°C before proceeding with PCR amplification. For PCR amplification, we used 1.25 µl IDT for Illumina Nextera DNA Unique Dual Indexes Set C, which contains two indexes premixed and 25 µl of Bioline MyFi Mix. This is in the place of the NEB Next HiFi PCR mix in your protocol. We performed the PCR as outlined. For PCR amplification, 15 cycles were used due to the reduced amount of material. The amplified library was purified using the Qiagen PCR purification MinElute kit. This was followed by a 1× volume AMPure XP bead-based clean-up according to manufacturer's guidelines. The resulting libraries were evaluated for quantity and quality using a Thermo Fisher Scientific Qubit 4.0 fluorometer with the Qubit dsDNA HS Assay Kit and an Advanced Analytical Fragment Analyzer System using a Fragment Analyzer NGS Fragment Kit, respectively.

The ATAC-Seq libraries were further quantified by qPCR using a Bioline JetSeq library Quantification Lo-ROX kit according to their guidelines. The libraries were pooled equimolarly and further cleaned using AMPure XP beads as described above. The library pool was then again assessed for quantity and quality using fluorometry and capillary electrophoresis, as described above.

The pool was loaded at 150 pM using an XP workflow into one lane of a NovaSeq 600 SP with NovaSeq XP 2-Lane Kit v1.5. The libraries were sequenced paired end on an Illumina NovaSeq 6000 sequencer using a NovaSeq 6000 SP Reagent Kit v1.5. An average of 56 million reads/library were obtained. The quality of the sequencing runs was assessed using Illumina Sequencing Analysis Viewer and all base call files were demultiplexed and converted into FASTQ files using Illumina bcl2fastq conversion software v2.20. ATAC-seq experiments were performed in collaboration with the Genomics Unit of the University of Bern.

All bioinformatics analyses were performed using bash scripts or R statistical software. A quality check of the samples was performed using FASTQC and reports were summarized using MultiQC (Andrew, 2010; Ewels et al., 2016). Adapters from the fastq files were trimmed using trimmomatic software (Bolger et al., 2014). Reads were aligned using bowtie2 (Langmead and Salzberg, 2012; Langmead et al., 2019) to GRCz11 danRer11 v102 assembly from Ensembl (Yates et al., 2020) with flags ‘--very-sensitive’. Paired-end reads were used for downstream analysis. The files were then converted to bam, downsampled to the lowest counts and indexed, and the mitochondrial chromosome was removed using samtools (Danecek et al., 2021; Li et al., 2009). Duplicates were removed using Picard tools (https://github.com/broadinstitute/picard). The samples were processed to select for unique reads using samtools (Danecek et al., 2021; Li et al., 2009). The peaks were identified using Genrich in ATACSeq mode and zebrafish genome size (Gaspar, 2018). The zebrafish genome size was estimated using faCount script from public utilities from UCSC (Kent et al., 2002).

To analyze the differential accessible regions, we used DiffBind using background and DeSeq2 normalization and a cutoff threshold of P<0.05. To annotate the peaks, we used ChIPSeeker (Yu et al., 2015). We used the GRCz11 danRer11 v102 assembly from Ensembl and transcription start site region as ±1 kb for annotation. The annotated genes were then converted to mouse orthologous genes using biomaRt and used for pathway analysis using clusterProfiler (Durinck et al., 2009; Yu et al., 2012). K-means clustering was performed using SeqMINER software using linear enrichment clustering approach with 10 clusters (Ye et al., 2011). The bigwig files to visualize the peaks were made using bamCoverage in deepTools2 (Ramírez et al., 2016). Interactive Genome Viewer was used to visualize the peaks (Robinson et al., 2011).

To identify the transcription factor-binding sites, we used the sequences from the differential accessible regions in Centrimo from MEME-suite. We used CIS-BP 2.0 Danio rerio Database to identify the potential zebrafish transcription factors (Bailey and Machanick, 2012; Weirauch et al., 2014). The detailed ATAC-seq pipeline used can be found at https://github.com/MercaderLabAnatomy/PUB_Marques_et_al_2021/commit/d40225b5751ab2089d9fff47f46f898fdcb1dfbf.

Immunofluorescence images were acquired using the Leica TCS SP8 DLS confocal microscopes. For image acquisition of whole-mount embryos, larvae were mounted in 1% low melting agarose in a MatTek Petri dish. Images were acquired with a 20× water immersion objective. Images were processed afterwards with Fiji software. Figure legends indicate whether a 3D projection is presented or a maximum intensity projection of a reduced number of stacks is shown. For 3D projections, images were first treated with a mean filter, with a radius of 2.0 pixels. Interpolation was also applied when rendering the 3D projections.

To assess the eGFP/mRFP ratio from in vivo confocal images, we applied a median filter (3 pixels radius) and measured line profiles from the SV 60 µm into the atrium in six sequential z-planes. The mean intensity along the line profile normalized by the maximum per fluorophore per embryo was calculated. Subsequently the ratio for each µm along the line profile was obtained.

For quantification of mean fluorescence intensity, a mean filter with a radius of 2.0 pixels was first applied to smoothen the images. Afterwards, we created a maximum intensity projection of all the stacks containing the heart. We then delimited the heart and applied an automatic OTSU threshold. Automatic threshold was evaluated independently for each image and, when necessary, minor adjustments were applied. Finally, mean fluorescence intensity was calculated.

Semi-quantification of signal intensity for whole-mount in situ hybridization was carried out using Fiji software. First the images were inverted and the region of interest (ROI) was defined and used for all images. For each image, mean signal was measured in six independent areas: three in the background and three in the stained area. Measurements were averaged and then background signal subtracted from the signal measured in the stained area. The fold change was calculated and GraphPad was used for statistical analysis.

Statistical analysis was carried out with GraphPad Prism 7. When data fitted normality parameters, i.e. passed either the D'Agostino-Pearson or the Shapiro-Wilk normality test, an unpaired t-test was used. If this was not the case, a Mann-Whitney non-parametric test was used to compare differences between conditions. If a statistically significant difference in the standard deviation between conditions was detected, the unpaired t-test with Welch's correction was applied. In case of multiple comparisons, a one-way ANOVA was applied, followed by Tukey's multi-comparison test. For each presented graph, the type of statistical test applied is stated in the figure legend.

We thank Anna Gliwa and Eduardo Diaz for fish husbandry at the University of Bern and Centro Nacional de Investigaciones Cardiovasculares, respectively. Microscopes supported by the Microscopy Imaging Center (MIC) at the University of Bern were used. We thank Stephan Müller from the FACS Lab of the University of Bern for help with FACS. We thank the Genomics Unit from CNIC for help with SMARTer-seq, and Pamela Nicholson and Cátia Coito from the Genomics Unit of the University of Bern, for help with ATAC-seq. We thank Didier Stainier for sharing the Tg(myl7:cdh2-tdTomato)^bns78 line and Julien Vermot for comments on the manuscript.