Transcriptome analysis of mouse and human sinoatrial node cells reveals a conserved genetic program

The sinoatrial node (SAN) is the dominant pacemaker of the heart, located at the junction of the right atrium (RA) and the superior caval vein near the crista terminalis. The SAN consists of pacemaker cells that are able to depolarize spontaneously, thereby initiating the cardiac cycle. The SAN is innervated by the autonomic nervous system, which enables modulation of the beat rate and output. Fibrosis, aging or loss of pacemaker cells may cause SAN dysfunction, resulting in bradycardia and sudden death (Adan and Crown, 2003; Dobrzynski et al., 2007). Ultimately, SAN dysfunction requires implantation of an electronic pacemaker. The generation of bona fide pacemaker cells by programming stem cells or reprogramming endogenous non-pacemaker cells has stimulated much interest, as biological pacemakers may overcome the shortcomings of the currently used electronic devices (Boink et al., 2015; Cingolani et al., 2017; Protze et al., 2017). The ability to create biological pacemakers may be improved by the molecular and developmental characterization of the SAN.

Compared with working cardiomyocytes, SAN cells contain poorly developed sarcoplasmic reticula and sarcomeres (Masson-Pevet et al., 1979; Op ‘t Hof, 1988). Furthermore, the electrical activity of SAN cells differs from working myocardium (WM) as a result of differences in expression of genes encoding ion-handling (e.g. channels, accessory units) and other proteins (Dobrzynski et al., 2007; Marionneau et al., 2005; van Kempen et al., 1991). For example, the center of the SAN conducts the impulse slowly and expresses low levels of Gja5 (Cx40), Gja1 (Cx43) and Scn5a (Na_v1.5), a set of genes that provide rapid conduction between atrial cardiomyocytes (Dobrzynski et al., 2007; Marionneau et al., 2005). Expression of Gjb6 (Cx30), Gjd3 (Cx30.2) and Gjc1 (Cx45) in the SAN has also been described (Gros et al., 2009; Kreuzberg et al., 2005). Compared with WM cells, SAN cells have a less-negative diastolic membrane potential, a slower action potential upstroke velocity and do not possess a stable resting membrane potential because of the absence of the inward rectifier K⁺ current (I_K1), which is carried by the K_ir2 channel family (Chandler et al., 2009). In addition, SAN cells exhibit the hyperpolarization-activated current, which is an inward current carried by K⁺ and Na⁺ ions (Brown et al., 1979), and is also referred to as ‘funny current’ (I_f). This current may contribute to the diastolic depolarization phase of SAN cell action potentials (DiFrancesco, 2006; Verkerk et al., 2009), and is carried by the hyperpolarization-activated cyclic-nucleotide gated (Hcn) gene family, of which Hcn4 is the most predominant isoform found in the SAN (Marionneau et al., 2005; Moosmang et al., 2001; Stieber et al., 2003). Owing to its virtual absence in the WM and its abundance in the SAN area, Hcn4 has become a potent marker for the mature SAN, although its expression in the mouse and human embryonic hearts, and during pathophysiological conditions, is much broader (Liang et al., 2015a; Mommersteeg et al., 2007; Nattel et al., 2008; Sizarov et al., 2011). Mutations in HCN4 in humans are associated with bradycardia (den Hoed et al., 2013; Milano et al., 2014; Verkerk and Wilders, 2015). Moreover, Hcn4-deficient mice die in utero due to severe bradycardia (Stieber et al., 2003). In addition to ion channels, pacemaker activity is also regulated by a tight coupling of sarcoplasmic reticulum (SR) Ca²⁺ cycling molecules with the electrogenic Na-Ca exchanger, also known as the ‘Ca²⁺ clock’ (Lakatta et al., 2010; Mangoni and Nargeot, 2008).

Previous studies have shown that core transcription factors are involved in SAN formation (reviewed by van Eif et al., 2018), including Tbx3 (Hoogaars et al., 2007; Mommersteeg et al., 2007; Wiese et al., 2009), Isl1 (Liang et al., 2015b; Tessadori et al., 2012), Shox2 (Blaschke et al., 2007; Hoffmann et al., 2013; Puskaric et al., 2010) and Tbx18 (Mommersteeg et al., 2010; Wiese et al., 2009). Moreover, the key cardiac transcription factor Nkx2-5 is selectively repressed in the developing SAN (Espinoza-Lewis et al., 2009; Mommersteeg et al., 2007; Ye et al., 2015). BMP and WNT signaling might contribute to SAN development (Bressan et al., 2013; Hashem et al., 2013; Puskaric et al., 2010). Lineage tracing reveals that the definitive SAN develops from Tbx3⁺ cells in the early heart tube (Mohan et al., 2018). Tbx3 has been shown to activate the pacemaker gene program and is able to impose pacemaker function on the atria (Bakker et al., 2012; Hoogaars et al., 2007). Furthermore, Tbx3 suppresses genes associated with atrial WM (Bakker et al., 2012; Hoogaars et al., 2007; Singh et al., 2012). Thus, Tbx3 is a potent regulator and marker of the SAN, and can be considered as a suitable target gene for functional analysis of pacemaker tissues.

To gain insight into the molecular mechanisms underlying SAN regulation and function, we set out to obtain the transcriptome of the SAN pacemaker cells and to measure functional parameters of the pacemaker cells throughout development. However, the SAN comprises a heterogeneous cell population, consisting of central pacemaker cells, peripheral pacemaker cells, fibroblasts, endothelial cells and connective tissue (Dobrzynski et al., 2007; Sanchez-Quintana et al., 2005). Because this heterogeneous cellular composition interferes with obtaining uncontaminated pacemaker cardiomyocytes for analysis, we generated a mouse model that expresses the fluorescent protein Venus under control of the Tbx3 locus in the Tbx3⁺ cells. Patch clamp and optical mapping experiments on isolated Venus⁺ and Venus⁻ right atrial cells at embryonic, fetal and adult stages were performed to study the electrophysiological properties of pacemaker cells during development, and to validate the mouse model. To determine the gene transcriptional profile of the SAN, we obtained the transcriptome of FACS-purified Venus⁺ SAN cells and Venus⁻ working cardiomyocytes. In addition, the transcriptome of fetal human SANs was assessed and compared with that of mouse, in order to identify conserved genetic networks driving the pacemaker gene program. We identified Smoc2 as a novel conserved highly specific SAN marker, and inactivated the gene to assess its function.

In order to visualize the cardiac conduction system, the reporter gene encoding the yellow fluorescent protein Venus was inserted into the start codon of Tbx3 (Fig. 1A and Fig. S1A). Whole-mount in situ hybridization revealed expression of Tbx3 and Venus in the snout, eye, ear, mammary glands, genital tubercle and limbs of E10.5 wild-type and Tbx3^+/Venus embryos, respectively (Fig. 1A). The pattern of fluorescent Venus signal was highly similar to that of endogenous Tbx3 (Fig. S1D). Because the Venus-coding sequence was inserted into the translation start site of Tbx3, expression of Tbx3 from the modified allele is expected to be disrupted. To validate inactivation, we generated homozygous mutants (Tbx3^Venus/Venus) by intercrossing heterozygous Tbx3^+/Venus mice. Tbx3^Venus/Venus embryos died before developmental stage E12.5 and exhibited abnormal forelimbs, a double outlet right ventricle and failed liver development (Fig. S1E), in accordance with the phenotype of Tbx3 homozygous mutants described previously (Bakker et al., 2008; Davenport et al., 2003; Frank et al., 2011; Lüdtke et al., 2009; Mesbah et al., 2008). Tbx3^+/Venus mice, in contrast, did not show any abnormalities and were viable and fertile.

To investigate whether Venus is expressed selectively in the Tbx3⁺ SAN, immunohistochemistry was performed on embryonic, fetal and adult hearts. Venus was co-labeled with antibodies labeling the SAN (Tbx3 and Hcn4), atrial working cardiomyocytes (Gja5) and cardiomyocytes in general (cTnI). We found that cTnI⁺/Hcn4⁺/Tbx3⁺ cells co-expressed Venus, which indicates SAN pacemaker cardiomyocytes are accurately identified by Venus expression (Fig. 1B and Fig. S2). No Venus expression was detected in the Cx40⁺/Hcn4⁻/Tbx3⁻ atrial working cardiomyocytes. Within the SAN region, Venus was not detected in cTnI⁻/Hcn4⁻/Tbx3⁻ non-cardiomyocytes. However, Venus expression was seen in the Tbx3⁺ cardiac ganglia, which is located at the dorsal side of the atria near the SAN. These results imply that Venus is expressed in a pattern identical to that of Tbx3 and that Tbx3^+/Venus mice are suitable to mark the SAN pacemaker cardiomyocytes.

To gain insight into the structure of the SAN, we performed a three-dimensional reconstruction of the Venus⁺ area and atria (Fig. 1C). Although the SAN is usually depicted as a comma shaped structure (Boyett et al., 2000; Wiese et al., 2009), the Venus⁺ area appears as a ring structure around the superior caval vein (SCV). The SAN ‘head’ is found at the cranial and right side of the SCV, and is connected with the ‘tail’ running from the sulcus terminalis to the crista terminalis. In addition, Venus⁺ myocardium forms a connection between the left and right atrium and is further visible at the left side of the SCV in the interatrial septum, connecting the SAN with the atrioventricular conduction axis. The existence of both left and right tracts – known as the anterior (left) and posterior (right) intermodal pathway – has been described in mouse, rat and human previously (Aoyama et al., 1993; Blom et al., 1999; James, 2001; Rentschler et al., 2001). To investigate whether SAN size is affected by reduced Tbx3 expression in Tbx3^+/Venus mice, we determined the shape and volume of Hcn4⁺ SAN regions in ND0 wild-type and Tbx3^+/Venus mice (both N=3). Comparison revealed that Tbx3 heterozygosity does not result in altered SAN shape or size (Fig. 1D).

To study whether the site of initial activation of the atria overlaps with the Venus-expressing (SAN) domain, we performed high-density optical mapping on Tbx3^+/Venus tissue preparations (n=14), containing the RA, crista terminalis and intercaval area. Mapping the endocardial (n=8) and epicardial (n=6) sides of the preparations revealed initial pacemaker activity in the center of the Venus⁺ domain (Fig. 2A). To test sympathetic activity, samples were exposed to 1 µM isoproterenol. Upon stimulation, beating rates increased from 395±27 beats per minute (bpm) to 542±23 bpm (P<0.05), indicating that the sensitivity of the SAN to adrenergic stimulation is maintained in Tbx3^+/Venus mice (Fig. S3A).

Patch-clamp experiments were performed on isolated cells to investigate the electrophysiological properties. APs of adult Venus⁺ cells and isolated wild-type SAN cells showed no differences in cycle length, RMP, AP duration and DDR, indicating that the heterozygous loss of Tbx3 did not affect the function of the SAN cells of Tbx3^+/Venus mice (Fig. S2B,C). Venus⁻ cells in the SAN region exhibited APs associated with WM. In contrast, Venus⁺ cells showed spontaneous activity and a diastolic depolarization phase (Fig. 2B). Moreover, Venus⁺ cells displayed a significantly less negative RMP, a lower APA and upstroke velocity, and shorter AP durations at 20% and 50% of repolarization than Venus⁻ cells (Fig. S3B). In addition, we detected an I_f current in Venus⁺ cells, whereas a Na⁺ current (I_Na; encoded by Scn5a) was not present (Fig. S3C). To study pacemaker cell characteristics during development, we compared isolated Venus⁺ cells from E12.5, E17.5 and adult hearts. At subsequent stages of development, SAN pacemaker cells showed increasing automaticity, as displayed by the decrease in cycle length and increase in DDR (Fig. 2C-D). Furthermore, RMP became more negative and APA increased. Finally, APs became shorter and the upstroke velocity increased (Fig. 2C,D).

We obtained transcriptional profiles of both late fetal (E17.5) SANs and working cardiomyocytes using RNA-sequencing. To obtain purified SAN and working cardiomyocytes, we crossed Tbx3^+/Venus with BAC336-Nppb-Kat mice, in which red fluorescent protein Katushka is expressed selectively in the cells of the working cardiomyocytes (Sergeeva et al., 2014). Double transgenic mice (Tbx3^+/Venus;BAC336-Nppb-Kat) have Venus⁺/Katushka⁻ SAN and Venus⁻/Katushka⁺ working cardiomyocytes (Fig. 3A). The SAN, atria and ventricles of double heterozygous E17.5 fetuses were isolated, dissociated and FACS-purified based on Venus⁺ or Katushka⁺ fluorescence (Fig. 3B). Double-positive Venus⁺/Katushka⁺ cells were not observed. After performing RNA-seq on three samples of ±7000 (Venus⁺) and ±250,000 (Katushka⁺) cells each, principal component analysis (PCA) on both tissue populations revealed the working cardiomyocyte transcriptomes to be similar, and to be substantially different from the SAN (Fig. S4A). On the other hand, larger variation was observed between the SAN samples. The small SAN size and low cell numbers, along with incomplete cell dissociation before FACS purification, are possible explanations for sample heterogeneity. Although distinction between single cells and duplicates is possible, the manual gate setting is arbitrary (Fig. S4B). We identified 16,974 genes to be expressed in one or both cell types [normalized reads per kilobase (RPK)>10 reads] (Fig. S4C; Table S1).

To verify the purity of the SAN and working cardiomyocyte cell populations, we compared the expression patterns of SAN and working cardiomyocyte markers Tbx3, Shox2, Isl1, Gja5, Nppa and Nppb between both cells types (Blaschke et al., 2007; Hoffmann et al., 2013; Hoogaars et al., 2007; Liang et al., 2015b; Sun et al., 2006; Vedantham et al., 2015). Read counts for Tbx3 and Shox2 were 60- and 63-fold enriched, respectively, in the Venus⁺ cells (Fig. S4C). Higher Gja5, Nppa and Nppb transcripts levels (4-, 15- and 100-fold enriched, respectively) were found in the Katushka⁺ population (Fig. S4C). Based on the near absence of reads of working cardiomyocyte markers and very high read counts of SAN markers in each individual SAN transcriptome, we conclude that the quality of the SAN samples is adequate (Fig. S4C). We performed cluster analysis on significantly enriched genes, showing two main groups separated by cell type (Fig. 3C). At a significance level of P<0.05 and log 2-fold change <−0.5 or >0.5, we identified 3109 SAN-enriched genes and 1655 working cardiomyocyte-enriched genes (Fig. 3D,E). Besides the known pacemaker cell markers, we found a subset of Ca²⁺ channels (Cacna2d2 and Cacna1g), putative Ca²⁺-binding proteins (Smoc2 and Vsnl1), ligands and effectors of the BMP-signaling pathway (Bmp2, Bmp4 and Smad9) and neuron-associated genes (Th, Cadps and Nrp2) to be SAN enriched (Fig. 3D). Functional annotation on three main clusters revealed that genes expressed in mitochondria (Coq3) and contributing to the contractile apparatus (Ttn and Myh7) and fast conduction (Gja1 and Gja5) were working cardiomyocyte enriched (Fig. 3C,F). However, cells in the SAN region expressed higher levels of genes associated with neuronal function (Chrnb4), NOTCH- and BMP-signaling pathways, and calcium handling (Cacna2d2, Cacna1g) (Fig. 3C,F). Genes functionally annotated as ‘extracellular region’ and ‘cell adhesion’ were present in clusters containing working cardiomyocyte and SAN markers, although the latter showed highest significance for these annotation terms (6.5E-19 and 5E-22, respectively). A list of E17.5 SAN- and working cardiomyocyte-enriched genes is provided in Table S1.

We used in situ hybridization to validate the expression patterns of novel highly SAN region-enriched genes (Alox8, Smoc2, Vsnl1, Gfra2, Chrnb4 and Slc9a3r2) (Fig. 4A). Arachidonate 8-lipoxygenase (Alox8) was found to be expressed predominantly in the SAN head of E17.5 fetuses, whereas SPARC-related modular calcium binding 2 (Smoc2) and visinin-like protein 1 (Vsnl1) marked the entire SAN. Expression of Vsnl1 has been detected previously in the developing heart at the venous pole region, although its function remains unknown (Ola et al., 2012). Smoc2 has been shown to be involved in teeth development (Bloch-Zupan et al., 2011) and can be considered as an early intestinal stem cell marker (Muñoz et al., 2012). Interestingly, these novel highly pacemaker-enriched genes have not been associated with SAN function before. GDNF family receptor α2 (Gfra2) and neuronal acetylcholine receptor subunit β4 (Chrnb4) are expressed in the cardiac ganglia, although Gfra2 also showed expression, to a lesser extent, in the SAN head and tail. Furthermore, we found that sodium-hydrogen exchange regulatory cofactor 2 (Slc9a3r2) marks the epithelium in the ventricles and the SAN, including the nodal artery, but to a lesser extend in the atria.

The transcriptome analysis indicated BMP signaling to be SAN enriched. Using in situ hybridization, we found Bmp2 to be expressed in the Tbx3⁺/Nppa⁻ SAN region (Fig. 4B,D). Moreover, immunohistochemistry revealed co-staining of phosphorylated Smad1/5/8 and Tbx3 in the SAN, indicating elevated BMP-signaling activity in the SAN (Fig. 4C). To investigate whether BMP signaling affects the expression of SAN and WM genes, we incubated dissected fetal RA and LA WM explants and HL-1 cells with recombinant BMP2 for 48 h and performed qPCR analysis (n=3) (Fig. 3E-G). We found Cacna2d2 expression to be significantly induced in RA and LA explants after exposure to BMP2 (Fig. 4F). Moreover, incubation of dissected SANs with Bmp inhibitor Dmh1 lowered Cacna2d2 levels (Fig. 4F). We also analyzed the expression levels of other SAN and WM markers using RNA sequencing (Fig. 4G). Cacna2d2 and Cacna1g levels were significantly elevated, and Bmp2 and Tbx3 was slightly but not significantly induced. Conversely, WM markers were slightly reduced, and Scn5a was significantly decreased after activation of BMP signaling (Fig. 4G). These observations implicate a role for BMP signaling in activating the pacemaker gene program and the downregulation of WM genes.

Comparing transcriptomes of the SAN derived by independent preparation methods will improve the accuracy of identifying transcriptional profiles of the SAN. Furthermore, analysis of the gene expression profile at different stages provides insight into the mechanisms underlying SAN development. Therefore, we compared our E17.5 samples with the SAN transcriptome published previously (Vedantham et al., 2015) (Table S1). In this study, SANs were identified based on GFP expression in Hcn4⁺ cells and isolated using laser capture microscopy (LCM) at different stages (E14.5, P4 and P14). Because the SAN region comprises a mixture of different cell types (cardiomyocytes, epicardial cells, fibroblasts, endothelial cells, etc.), we expected that these samples would contain non-pacemaker cells within the LCM-isolated region. The E14.5 sample in addition could contain embryonic atrial cells, which still express some Hcn4 at this stage (Liang et al., 2013; Sizarov et al., 2011). Genes that were significantly SAN enriched (P<0.05) at one or more stages (E17.5 FACS, E14.5 LCM, P4 LCM and P14 LCM), which yielded 6472 genes, were clustered. This revealed that genes in clusters associated with platelet activation, fibroblast proliferation and neurons are expressed at higher levels in the E17.5 FACS-purified Venus⁺ SAN population compared with samples isolated by LCM (Fig. 5A). At all stages, genes involved in cardiac muscle contraction, respiratory chain function and myofibril assembly are more abundantly expressed in working cardiomyocytes. Effectors of the BMP-signaling pathway (Bmp2 and Bmp4) are SAN enriched in E14.5, E17.5, P4 and P14 mice. The SAN transcriptome is annotated with ‘negative regulation of myoblast differentiation’, whereas working cardiomyocytes express genes associated with positive regulation of cardiac muscle cell proliferation (Fig. 5A). To study the consistency between FACS- and LCM-derived samples in more detail, we assessed the relative expression levels of working myocardial, fibroblast, endothelium and neuronal cell markers in the populations (Fig. S5A-D). Normalized expression levels were corrected for normalized RPK values of a subset of known and validated pacemaker-specific genes (geomean of RPK of Isl1, Tbx3, Shox2, Smoc2, Bmp2 and Bmp4). Although the E14.5 SAN samples (LCM) contained lower levels of transcripts of fibroblast, endothelium and neuronal genes, the working cardiomyocyte marker Nppa is present at higher levels compared with other stages, indicating a lower degree of discrimination between SAN and working cardiomyocyte cells in these samples. Neuronal genes are more abundantly expressed in the E17.5 FACS Venus⁺ samples compared with LCM, which indicates the presence of more neuronal cells. As Tbx3 is expressed in cardiac ganglia (Horsthuis et al., 2009) (Fig. 1B), these cells may have been co-isolated in the Tbx3^+/Venus mouse model.

Because the FACS-purified and LCM-isolated samples differ in cell-type composition, we hypothesize that transcripts consistently enriched in both LCM and FACS SAN samples are more likely to be derived from pacemaker cardiomyocytes. We compared expression profiles of E14.5 LCM and E17.5 FACS samples, and found 2249 genes to be significantly higher expressed in the SAN in both samples (>2 fold expression) (Fig. 5B). Functional annotation revealed that genes involved in Ca²⁺ handling, NOTCH- and BMP-signaling pathway, and extracellular matrix organization were SAN enriched (Fig. 5C). In addition, we found a significant number of genes (>90) associated with neuronal expression. These genes may be derived from neuronal cells (cardiac ganglia) or could be expressed by the pacemaker cells themselves (Ebert and Thompson, 2001; Gorza et al., 1988; Horsthuis et al., 2009). To investigate whether the embryonic SAN gene program is maintained at adult stages, we dissected adult SANs and RA (n=2), derived from Tbx3^+/Venus;BAC336-Nppb-Kat mice, and found 13,500 genes with RPK-values of at least 10 in the SAN or RA WM tissues (Fig. S6A). Log2-fold changes from dissected adult samples (SAN/WM) were compared with E14.5 LCM (Fig. S6B) and yielded 379 genes, SAN enriched at both stages (L2FC<−1). Tbx3, Isl1, Shox2, Hcn4, Smoc2, Vsnl1 and BMP-signaling components were more abundantly expressed in the adult SAN, indicating that these genes robustly mark the SAN throughout all stages of life. Interestingly, Cntn2 was highly SAN enriched at the adult stage in contrast to the E14.5 stage, at which no expression was detected. We further found Cacna1d to be SAN enriched at both stages, and Cacna1b was more abundantly expressed at the adult stage. However, whether Cacna1b levels are derived from neuronal innervation rather than pacemaker cells remains unknown.

In order to identify transcriptional regulatory components possibly involved in SAN development or function, we performed cluster analysis on transcription factors (TFs) of E17.5 FACS, E14.5, P4 and P14 LCM samples that were significantly SAN-enriched at the 1 stage and over (Fig. S7A). For the SAN-enriched TFs, we found (among others) functional annotations referring to the BMP- and NOTCH-signaling pathways (Smads), and to cell fate commitment (Hoxa2). TFs that were enriched in all datasets (L2FC<0) are listed in Table S1. Besides known TFs (Tbx3, Isl1 and Shox2), SAN cells also express Hox-family members (Hoxa1, Hoxa2, Hoxa3, Hoxa4 and Hoxc4), BMP-signaling components (Smad1, Smad6 and Smad9) and Sox genes (Sox8 and Sox10).

Pacemaker cell-specific electrophysiological characteristics are caused by a subset of ion channels in the SAN, the expression of which is different during development, reflecting the alterations in APs (Fig. 2C,D). To characterize the ion channel distribution in the SAN and WM, we compared expression profiles of ion channel genes in E17.5 FACS-purified Venus⁺ and Katushka⁺ cell populations (Fig. S7B). SAN- and chamber-enriched ion channels are listed in Table S1. We identified 128 SAN-enriched ion channels (L2FC<0) (Fig. S7B). In contrast to a previous study (Marionneau et al., 2005), we found that Kcna5 (an ultrarapid delayed rectifier potassium channel), which is activated during the early repolarization phase, is SAN enriched rather than chamber enriched. As expected, we found Hcn4 and also Hcn1 to be more abundantly expressed in the SAN. In order to explain the differences in APs of pacemaker cells (Dobrzynski et al., 2007; Mangoni and Nargeot, 2008) during development, we determined the average expression levels of ion channels in E12.5, E14.5, P4 and P14 LCM samples. Expression of Atp2a2 [which encodes the sarco-endoplasmic reticulum (SR) Ca-ATPase SERCA2A) and Pln increased after birth. Expression of Ryr2 and Ryr3, which are responsible for sarco-endoplasmic reticulum Ca²⁺ release, was slightly increased. Expression of Hcn1 and Hcn4 (underlying the I_f current) also increased after birth (Fig. S7B), correlating with the developmental increase of phase 4 depolarization (DDR; Fig. 2). Expression of Scn5a, which underlies the main cardiac sodium current, was not altered, whereas Scn10a levels declined. Cacna2d2 and Cacna1c levels were not notably altered, but we found Cacna1g and Cacnb2 expression to be increased during development. Cacna1g (a T-type Ca²⁺ channel), in collaboration with I_f currents, is responsible for phase 4 depolarization towards the AP threshold. Cacnb2 (a L-type calcium channel) is responsible for the AP upstroke velocity in pacemaker cells, which was increased during development (Fig. 2C,D). It is also responsible for an important inward current during the AP; thus, the observed developmental AP shortening indicates upregulation of repolarizing K⁺ currents. Indeed, expression of Kcnd3, Kcnh2 and Kcnq1 (which underlie transient and rapid and slow delayed rectifier K⁺ currents) increases during development. We also found an increased expression of Kcnj3, Kcnj5, Kcnj11, Kcnj12 and Kcnj14, which are responsible for various inward rectifier K⁺ currents; this explains the developmental hyperpolarization of the RMP, despite the increased expression of Hcn4, Hcn1 and Cacna1g (Fig. 2).

To identify the transcriptome of human SANs, we isolated 9- and 12-week-old (both n=2) human fetal SAN regions and right atrial tissue by microdissection. As in mice, PCA analysis revealed that human pacemaker cells show a distinct phenotype compared with microdissected non-SAN right atrial tissue (Fig. S8A). We compared the log2-fold changes (SAN/RA non-SAN) of 9- and 12-week-old fetuses and found that 1076 and 698 genes, respectively, were SAN region and RA enriched (L2FC<−0.4) (Fig. 6A). ISL1, SHOX2, TBX3, VSNL1 and SMOC2 were enriched in the human SAN regions, as well as genes belonging to signaling pathways (BMP, NOTCH), anterior/posterior pattern specification genes, and genes associated with neuronal development and extracellular matrix organization (Fig. 6B). Moreover, expression of genes associated with the contractile apparatus and fast conduction was lower in the SAN area than in right atrial tissue, in accordance with the expression profile in mice. To further compare the SAN region transcriptomes of human and mice, we compared the human and E17.5 mouse transcriptional profiles and identified 855 genes to be SAN region-enriched in both datasets (Fig. 6C). Functional annotation analysis indicates that BMP- and NOTCH-signaling pathway components, as well as TFs ISL1, SHOX2 and TBX3, and HOX-family members in the SAN area are conserved between human and mouse. The potential role of HOX-family genes is further underlined by the association with retinoic acid response (Fig. 6D). We also compared expression profiles of E14.5 with human SAN regions, yielding 784 enriched genes, and found similar functional annotations (Fig. S8B,C). We found 89 genes to be SAN region enriched at all prenatal stages (mouse E14.5, E17.5 and human 9w+12w; Table S1). We validated the expression levels of SAN region-enriched genes in human RA and SAN samples by qPCR (N=4) (Fig. 6F) and confirmed enrichment of the core program ISL1, SHOX2 and TBX3 genes in the SAN as well as members of the BMP family (BMP4 and SMAD9), the NOTCH pathway (NOTCH3 and DLL1) and calcium-binding protein-encoding genes (SMOC2 and VSNL1). We performed a gene set enrichment and network analysis (Wang et al., 2017) on mouse E17.5/human fetal SAN region-enriched transcripts (Figs S9, S10 and Table S1). Regulation of development, neural differentiation, extracellular matrix, growth factor binding, calcium and ion channel activity were among the most significantly enriched biological processes and molecular functions.

Because Smoc2 levels are enriched in both murine and human SANs, and its function in the heart unknown, we studied whether inactivation of Smoc2 affects SAN function. Using CRISPR/Cas9, a 22 bp deletion in exon 3 was made to induce a frameshift mutation. In heterozygous mice, expression levels were reduced by 50% compared with wild types (Fig. S11A). In homozygous mice, only low mRNA levels could be detected, probably caused by nonsense-mediated decay (Frischmeyer and Dietz, 1999). Smoc2^−/− mice were found to be viable and fertile, and appeared to be unaffected. We investigated Smoc2^−/− mutants for alterations in cardiac electrophysiology. Both in vivo and ex vivo measurements revealed that Smoc2 depletion did not significantly affect heart rate (in vivo, both wild type and Smoc2^−/−, n=15; ex vivo, wild type, n=7; Smoc2^−/−, n=10) (Fig. S11B). Smoc2^−/− mice exhibit a trend towards increased RR intervals (Fig. S11C) or heart rate variability during ex vivo conditions (Fig. S11D), although these differences were not significant. These results indicate that the effect of Smoc2 inactivation on cardiac function is minimal.

We show that fluorescent reporter Venus targeted to the Tbx3 locus recapitulates the expression pattern of Tbx3, and can be used to identify and analyze pacemaker cardiomyocytes. Optical mapping, patch clamp analyses and morphometric analysis indicated that the electrophysiological properties and size of Venus⁺ and wild-type SANs were similar, indicating the presence of the reporter in Tbx3 causing heterozygosity is of minor influence. However, the innervating autonomic nerves also express Tbx3. Therefore, Tbx3 heterozygosity may affect neuronal function and, for example, heart rates during stress conditions. Owing to the small size of the SAN and its heterogeneous tissue composition, studying the molecular biology of pacemaker cells has been challenging. Endogenous Venus expression facilitates improved identification, isolation and purification of pacemaker cells. In recent years, the generation of bona fide pacemaker cells has gained interest. Generation of a biological pacemaker will benefit from transcriptome analysis of endogenous SANs. In this study, we used Tbx3^+/Venus;BAC336-Nppb-Kat transgenic mice to provide transcriptomes of E17.5 and mature FACS-purified Venus⁺ SAN and Katushka⁺ working cardiomyocytes. Functional annotation clustering was performed to characterize the cell populations. We found expression of genes associated with anterior/posterior pattern specification, extracellular matrix organization, neuronal development and BMP- and NOTCH-signaling pathways to be enriched in the SAN. We identified novel markers (Smoc2, Vsnl1 and Alox8) not previously associated with SAN development and function, and confirmed their expression pattern with in situ hybridization. Inactivation of Smoc2 by CRISPR/Cas9 genome editing did not result in altered cardiac function. However, no studies were performed during stress conditions, which could elicit an altered phenotype.

We found that the cluster containing the most highly SAN-enriched genes (Tbx3, Isl1, Shox2, etc.) was functionally associated with ‘extracellular region’ and ‘cell adhesion’. Interestingly, functional clustering in a previous study revealed that genes associated with these processes were downregulated in the SAN of E12.5 murine embryos when Isl1 was inactivated (Liang et al., 2015b). This indicates that genes involved in extracellular matrix organization and cell adhesion are regulated by the transcription factor Isl1 in the SAN and may contribute to its function. Moreover, ectopic expression of Tbx2 or Tbx3 in the developing atrial myocardium stimulates BMP and TGFβ signaling, and induces extracellular matrix formation around the Tbx2/3-expressing cardiomyocytes (Singh et al., 2012). This suggests that T-box factors in the pacemaker cells may function likewise. The extracellular matrix consists of collagen fibers and forms a regularly patterned framework in and around the SAN that limits intercellular contact and limits the conduction velocity (Boyett et al., 2000; Matsuyama et al., 2004; Sanchez-Quintana et al., 2005). Diseases in extracellular matrix organization are associated with cardiac arrhythmias (Adan and Crown, 2003; Thery et al., 1977).

In contrast to the SAN, WM cells highly expressed genes involved in muscle contraction, fast conduction and mitochondrial function. Indeed, the SAN contains fewer mitochondria and shows less contractile properties (Boyett et al., 2000; Dobrzynski et al., 2007; Virágh and Challice, 1980). We detected expression levels of Gja1, Gja5 and Scn5a in SAN samples that are not expressed in the nodal center. These transcripts may be derived from tissue with higher conduction velocity in the SAN periphery, reflecting its heterogeneity. Therefore, Scn5a⁺/Gja5⁺/Hcn4⁺/Nppa⁻ peripheral cells (Dobrzynski et al., 2005; Lei et al., 2004) may co-express (low) levels of Tbx3. Furthermore, Scn5a is expressed in the intercaval region and in larger cells in the SAN rather than in smaller cells (Maier et al., 2003). Single cell sequencing analysis may reveal the extent of heterogeneity of the pacemaker cells, all of which seem to express Tbx3 and Shox2. The fraction of Hcn4⁺ cells that co-express Isl1 is reduced after birth, indicating that Isl1 expression in the SAN becomes mosaic (Liang et al., 2015b).

We investigated whether the fetal E17.5 murine FACS-purified SAN transcriptome is comparable with LCM-derived transcriptional profiles of earlier stages of development and adulthood. We found that Isl1, Shox2, Tbx3, Hcn4, Smoc2 and Vsnl1, as well as the ion channels Cacna1g, Cacna2d2 and Chrnb4, are SAN enriched from stage E14.5 to adulthood. Interestingly, we found high levels of Cntn2 in the adult SAN samples, whereas no expression was detected at E14.5. Cntn2 was previously indicated to be expressed in cardiac conduction system components before birth (Pallante et al., 2010; Qiao et al., 2017). The expression profiling data implies that Cntn2 expression in the SAN is selective and starts after birth.

Contamination rates, i.e. the presence of non-SAN cardiomyocytes in the samples, may depend on the techniques used for cell isolation and purification. Fibroblasts, epicardial cells, endothelial cells, etc. are positioned in close proximity to pacemaker cells and might be separated more precisely by FACS purification than LCM. However, relatively high levels of fibroblast and endothelial genes were detected in P14 LCM and E17.5 FACS samples (Fig. S5). This could be explained by the gateway setting procedure prior to FACS purification, which is arbitrary (Fig. S4B). At E14.5, discrimination between SAN and right atrial cardiomyocytes was less evident due to higher Nppa expression levels. Both LCM-isolated and FACS-purified cell populations show abundant expression of neuronal genes, although this was more prominent in FACS-purified samples. The SAN is innervated by autonomic nervous system cells that play important roles in regulating cardiac function. As Tbx3 was detected in the cardiac ganglia (Fig. 1B) (Horsthuis et al., 2009), processing microdissected and FACS-purified Venus⁺ SAN cells is more susceptible to neuronal contamination. However, neuronal genes are known to be expressed in cells of the cardiac conduction system. Neurofilament was detected in the SANs of adult mice (Wen and Li, 2015) and in rabbits (Gorza and Vitadello, 1989; Verheijck et al., 1998). Neuronal gene expression was abundant in the atrioventricular canal (precursor of the atrioventricular conduction system) of E10.5 murine embryos, long before innervation occurs here (Horsthuis et al., 2009) and in the SAN and bundle branches of rats from E12.5 onwards (Nakagawa et al., 1993). Moreover, immunohistochemistry on rabbit heart sections with neurofilament marker iC8 revealed staining near the atrioventricular junction (Gorza et al., 1988). Furthermore, Pnmt that converts noradrenalin into adrenalin, is expressed in the SAN of E11.5 rats, days before cardiac innervation occurs (Ebert and Thompson, 2001). This implies that nervous tissue and developing pacemaker myocardium express an overlapping gene program and that the heart is capable of synthesizing and storing catecholamines itself.

Knowledge regarding the molecular biology of pacemaker cells, which is essential for the generation of biological pacemakers, relies on functional studies in mice and cell lines. Insight into the transcriptome of human SANs will improve the robustness of pacemaker cell therapies for application in humans. Owing to the small size of the SAN and lack of clear morphological markers, isolation of fetal human SAN cells is arbitrary. Nevertheless, induced pluripotent stem cells (iPSC) used for generating biological pacemakers are immature and therefore transcriptome analysis of fetal SAN cells is relevant. We performed transcriptome analysis on microdissected human fetal SAN region tissue of 9 and 12 weeks of gestation, and non-SAN right atrial tissue. We observed that a gene program is conserved between human and murine SAN regions, including transcription factors (ISL1, SHOX2, TBX3, Bmp and Notch effectors and HOXA-family genes) and calcium-binding proteins (SMOC2 and VSNL1). In addition, we found genes associated with neuronal development to be SAN region enriched in both species. Because microdissection of human fetal pacemaker cells is challenging, a relatively high fraction of atrial and non-pacemaker cells was found in the SAN samples, explaining the prominent differences in gene transcriptomes between 9- and 12-week-old SAN regions (Fig. S7), and the fairly low number of genes consistently expressed at significant higher levels in both SAN area preparations.

We found Notch and Bmp effectors to be SAN region-enriched in human and mice. Tbx5, Shox2 and Bmp4 were found to play a role in inflow tract development of the mouse embryos (Puskaric et al., 2010) Isl1 deficiency resulted in reduction of Bmp4 expression in the SAN region (Liang et al., 2015b; Vedantham et al., 2015). The role of Bmp4 in SAN development remains elusive. Our data indicate specific expression of Bmp2/4 ligand and elevated BMP signaling in the developing SAN of mouse and human. Canonical WNT signaling might be involved in proliferation of the Tbx18⁺ mesenchymal progenitor cell population at the venous pole. The dorsal part of the sinus horns of Ctnnb1-deficient mice was not myocardialized (Norden et al., 2011). However, SAN formation was unaffected (Norden et al., 2011). Canonical WNT signaling was shown to be active in the SAN and AVC region in E10.5 murine embryos (Gillers et al., 2015). In mouse models, a loss-of-function mutation in WNT signaling is associated with septum and tricuspid valve defects. Abnormal chamber development and a prolonged PR interval and QRS complex have also been described (Gillers et al., 2015).

Transcription factors belonging to the Hox family were identified as SAN region enriched in both murine and human SANs. It has been suggested that initial Hox gene expression is activated by retinoic acid (RA) signaling, as anterior Hox genes exhibit abnormal expression patterns in embryos with impaired RA synthesis (Raldh2^−/− mutant) (Niederreither and Dollé, 2008; Waxman and Yelon, 2009). In zebrafish, Hoxb5b acts downstream of RA signaling and restricts the number of cardiac progenitor cells (Waxman et al., 2008). Moreover, expression patterns of Hoxa1, Hoxa3 and Hoxb1 show overlap with those of Raldh2 in E7.25 murine embryos (Bertrand et al., 2011). Mice lacking Raldh2, encoding the enzyme responsible for synthesizing the majority of embryonic RA, display an impaired expansion of second heart field progenitor cells and morphogenetic defects at the arterial and venous poles (Ryckebusch et al., 2008). The Hox family genes Hoxa1, Hoxa3 and Hoxb1 may be involved in cardiac second heart field development (Bertrand et al., 2011). First, genetic lineage tracing analysis has shown that Hoxb1⁺ and Hoxa1⁺ cells give rise to the inferior wall of the distal OFT and the atrial myocardium (Bertrand et al., 2011). Interestingly, development of the atria and sinus venosus is severely impaired in Raldh2 mutants (Niederreither et al., 2001). Whether this phenotype is caused by depletion of Hox gene expression has not been fully established.

Animal care was in accordance with national and institutional guidelines. Human fetuses (9 and 12 weeks) were obtained from electively terminated pregnancies, anonymously donated to research after informed written consent from donors was obtained in concordance with French legislation (PFS14-011) and with prior approval of the protocol (to S.Z.) from the ‘Agence de la biomédecine’.

Generation of BAC336-Nppb-Kat mice has been described previously (Sergeeva et al., 2014). A cosmid with Tbx3, isolated from the 129/Ola cosmid genomic library obtained from the Resourcenzentrum (RZPD) in Berlin, was kindly provided by Dr Andreas Kispert (Institut für Molekularbiologie, Medizinische Hochschule Hannover, Germany). Homologous DNA sequences (6.1 kb of upstream sequence and 1.9 kb of downstream sequence) were ligated to a Venus-polyA-Frt-flanked PGK-neo cassette derived from pKOII (Bardeesy et al., 2002) to generate a Tbx3-targeting construct (Fig. S1A) in which the first three codons of the Tbx3-coding region were replaced by the Venus-pA cassette (Nagai et al., 2002). The linearized targeting construct was electroporated into E141B10 embryonic stem (ES) cells to generate targeted cell lines. A diphtheria toxin A cassette was used to positively select for homologous recombinants. Chimeras were generated by injection of targeted ES cells into C57Bl6 host blastocysts. Germline transmission of the targeted allele was obtained by mating with FVB females. Subsequently, Tbx3-VenusNEO mice were crossed with FlpE mice (Rodríguez et al., 2000) to remove the PGK-neo cassette. Removal of the Neo cassette was assessed by Southern blotting and PCR (Fig. S1B,C). Progeny were screened by PCR for the presence of the Tbx3-Venus allele using the following primers: fw1 (AGCGGAGCCAAGCCAGCA), rv1 (CCTTGGCCTCCAGGTGCAC) and rv2 (TTGATGCCGTTCTTCTGCTTGT). The Tbx3-Venus allele has been maintained on the FVB genetic background.

DNA target sites and corresponding oligos were identified using ZiFit Targeter software. Oligos were annealed [95°C for 5 min and 25°C for 1 min (ramp down 0.1°C/s)] using 1× T4 ligation buffer (Invitrogen) and ligated into BsaI-digested and gel purified sgRNA expression vector pDR274, using T4 DNA ligase (Invitrogen). DNA was isolated with the JETSTAR DNA isolation kit (ITK, 200050). Expression vectors for sgRNAs and Cas9 (3 µg each) were linearized for 4 h at 37°C with DraI (NEB, R0129S) and PmeI (NEB, R0560S), respectively, followed by phenol-chloroform purification. Transcription of sgRNAs and Cas9 was performed using the MEGAshort T7 kit (Life Technologies, AM1354) and mMessage mMachine T7 Ultra kit (Life Technologies, AM1345), respectively, followed by purification with the MEGAclear kit (Life Technologies, AM1908), according to the manufacturer's protocol. For the injections, we used 5 ng/µl sgRNA and 10 ng/µl Cas9 mRNA.

Single cells were isolated from the SAN region and right atria of E12.5, E17.5 or adult hearts by an enzymatic dissociation procedure modified from Verkerk et al. (2009). In brief, the SAN region and atria were excised, cut into small strips (0.3-0.5×1 mm) and stored in a cold (6°C) modified Tyrode's solution containing (in mM): NaCl 140, KCl 5.4, CaCl₂ 1.8, MgCl₂ 1.0, glucose 5.5 and HEPES 5.0 (pH 7.4). Next, the strips were placed in nominally Ca²⁺-free Tyrode's solution (20°C), i.e. modified Tyrode's solution with 10 µM CaCl₂, which was refreshed twice. Then, the strips were incubated for 11-13 min in nominally Ca²⁺-free Tyrode's solution (37°C) to which liberase IV (0.25-0.29 U/ml, Roche) and elastase (2.4-0.7 U/ml, Serva) were added. During the incubation period, the strips were triturated through a pipette (tip diameter 2.0 mm). The dissociation was stopped by transferring the strips into a modified Kraft-Brühe (KB) solution (37°C), that was refreshed three times (at 20°C). The modified KB solution contained (in mM): KCl 85, K₂HPO₄ 30, MgSO₄ 5.0, glucose 20, pyruvic acid 5.0, creatine 5.0, taurine 30, β-hydroxybutyric acid 5.0, succinic acid 5.0, BSA 1% and Na₂ATP 2.0 (pH adjusted to 6.9 using KOH). Subsequently, the strips were stored for 30 min and then triturated (pipette tip diameter: 1.2 mm) in modified KB solution (20°C) for 2 min to obtain single cells. Single cells were stored at room temperature for at least 45 min in modified KB solution before they were put into a recording chamber on the stage of an inverted microscope, and superfused with modified Tyrode's solution (37°C). Individual Venus⁺, Venus⁻ or wild-type cells were isolated, and individual spindle and elongated spindle-like cells displaying regular contractions for SAN experiments were selected.

Membrane potentials and currents were recorded using the amphotericin-perforated patch-clamp technique using an Axopatch 200B amplifier (Molecular Devices). Signals were low-pass filtered at 10 kHz cut-off frequency, and digitized at 25 kHz. Data acquisition and analysis were accomplished using custom software. Pipettes (borosilicate glass; resistance 2–4 MΩ) were filled with a solution containing (in mM): K-gluconate 125, KCl 20, NaCl 5, amphotericin B 0.22 and HEPES 10 (pH adjusted to 7.2 using KOH). Action potentials (APs) were measured in both spontaneous active and quiescent cells. In quiescent cells, APs were elicited at 2 Hz by 3 ms current pulses, 1.5× threshold current pulses through the patch pipette. We analyzed cycle length, maximal diastolic potential (MDP), diastolic depolarization rate (DDR, measured over the 50 ms time interval starting at MDP+1 mV), maximal upstroke velocity (V_max), maximal AP amplitude (APA) and AP duration at 20, 50 and 90% repolarization (APD₂₀, APD₅₀ and APD₉₀). Parameters from 10 consecutive APs were averaged.

The mice were killed by cervical dislocation, after which the heart was excised and placed in a solution (30°C) containing (in mM): NaCl 128, KCl 4.7, CaCl₂ 1.45, MgCl₂ 0.6, NaHCO₃ 27, NaH₂PO₄ 0.4 and glucose 11 (pH was maintained at 7.4 by equilibration with a mixture of 95% O₂ and 5% CO₂). Fat and non-cardiac tissue was removed from the venous pole of the heart. Then the right atrium, crista terminalis and intercaval area were dissected free, whereafter the preparation was pinned down. Subsequently, the isolated preparations were incubated in 10 ml Tyrode's solution containing 15 μM Di-4 ANEPPS, superfused and placed in the optical mapping setup. Excitation light was provided by a 5 W power LED (filtered 510±20 nm). Fluorescence (filtered>610 nm) was transmitted through a tandem lens system on CMOS sensor (100×100 elements, sampling rate 5 kHz, MICAM Ultima). Optical action potentials were analyzed using custom-made software.

Data are mean±s.e.m. Two groups were compared using an unpaired t-test and more than two groups were compared using one-way ANOVA followed by a Student-Newman-Keuls method post-hoc test. P<0.05 indicates statistical significance.

Animals were anaesthetized with 5% isoflurane (Pharmachemie) and maintained under 1.5-2.0%. Electrodes were placed in the right (R) and left (L) armpit, and in the left groin (F). An electrocardiogram (ECG) was recorded (PowerLab 26T; AD-Instruments) for a period of 5 min. ECG parameters were determined in Lead II (L-R) based on the last 30 s of the recording. To record ECGs ex vivo, the adult mice were killed by cervical dislocation. The heart was rapidly excised, cannulated and mounted on a Langendorff perfusion set-up, as described previously (Boukens et al., 2013). Hearts from neonatal mice (ND0-1) were isolated and superfused with HEPES-buffered Tyrode's solution containing (in mM): 140 NaCl, 5.4 KCl, 1.8 CaCl₂, 1.0 MgCl₂, 5.5 glucose and 5.0 HEPES at a temperature of 36±0.2°C; pH was adjusted to 7.4 using NaOH. During perfusion, the heart was submerged and electrodes were placed at the right (R) and left (L) side of the base of the heart and at the left side of the apex (F) at a 5 mm distance. Lead II was used to determine ECG parameters.

E17.5 SANs and atrial WM tissues were dissected and placed in M199 medium (Sigma-Aldrich, M4530) supplemented with 2% fetal calf serum (Invitrogen, 10270-106), 1% penicillin/streptomycin (Thermo-Fisher, 15140-122) and incubated at 37°C. HL-1 cells were cultured at 37°C in Claycomb medium (Sigma-Aldrich, 51800C) containing 1% HL-1 defined FBS (Lonza, LO77227), 1% penicillin/streptomycin (Thermo-Fisher, 15140-122), 0.2 mM norepinephrine (Sigma-Aldrich, A0937-1G) (dissolved in 30 mM ascorbic acid solution) and 1% glutamax (Thermo-Fisher, 35050-061). HL-1 cells were seeded on plates coated with 0.02% gelatin and 0.5% fibronectin (Sigma-Aldrich, F-1141). HL-1 cells were recently authenticated and tested for contamination. Recombinant-Human BMP2 was added (R&D Systems, 355-BM) to atrial myocardial explants and HL-1 cells to a final concentration of 200 ng/ml. SAN samples were incubated with 3 µM Dmh1 (Sigma-Aldrich, D8946-5MG). All tissues were incubated for 48 h and medium was refreshed after 24 h.

Dissected tissue was collected in cold PBS (6°C) and dissociated using Trypsin/EDTA (0.25%/1 mM) to obtain a single cell suspension. Cells were passed through a cell strainer to exclude debris and cell clumps. Fluorescence-activated cell sorting (FACS) was performed on a Sony cell sorter (SH800Z, Sony Biotechnology). The human SAN and working atrial myocardium were recognized by their anatomic landmarks. Under the microscope, the right atrium was opened to expose the crista terminalis, the intercaval area and the interatrial septum. A thin strip of SAN tissue, limited by crista terminalis, atrial septum and orifices of the venae cavae was cut from the right atrium.

RNA from isolated SANs, atrial WM cells and cultured explants was obtained using the Nucleospin RNA XS isolation kit (Bioke, MN 740902.50), according to manufacturer's protocol. The quality of total RNA was assessed by micro-electrophoresis on acrylamide gel (Agilent 2100 Bioanalyser). For real-time RT-PCR Superscript II (Invitrogen) and OligodT primers were used to generate cDNA templates. Expression levels were assessed with quantitative RT-PCR using the Lightcycler Real-Time PCR system (Roche Diagnostics). Hprt was used as reference gene. For RNA-seq, the RNA template was processed with the Ovation V2 RNA-seq system (Nugen, 7102-32), according to the manufacturer's protocol. Libraries were generated using the Ovation Ultra-low system V2 1-96 (Nugen, 0344NB-A01). Amplified cDNA was bead purified (AmpureXP, Beckman-Coulter). Samples were sequenced on a Hiseq 2500 sequencing system (Illumina) with 50 bp single-end reads. RNA-seq data have been deposited in GEP under accession number available under GSE125932.

Hearts were fixed overnight in 4% formaldehyde, embedded in paraplast and sectioned at 7 μm. Embedding media were removed from sections and, by boiling for 5 min in a high pressure cooking pan in Antigen Unmasking Solution (Vector H3300), the epitopes in the tissue were unmasked. Tissue sections were incubated with primary antibodies overnight in Tris-NaCl buffer with blocking powder. Primary antibodies used were mouse-anti-troponin 1 (Millipore, MAB1691, 1:400), rabbit-anti-cardiac troponin 1 (Hytest LTD, 4T21_2, 1:200), mouse-anti-Cx40 (USBiological, C7856, 1:200), goat-anti-GFP (Abcam, 839963, 1:100), rabbit-anti-GFP (SZ E2110, 1:100), rabbit-anti-Hcn4 (Chemicon, NG164345, 1:200), goat-anti-Tbx3 (Santa Cruz, B1006, 1:200) and rabbit-anti-phospho-SMAD1/-5/-8 (Cell Signaling, 9511, 1:200). Secondary antibodies used were Alexa 680 (1:250; Molecular Probes, A21084), Alexa 568 (1:250; Molecular Probes, A10042) and Alexa 488 (1:250; Molecular Probes, A11017) with the appropriate epitope to visualize the primary antibody. For detection of goat anti-Tbx3 we used biotinylated donkey anti-goat antibody (1:250). For detection of biotinylated antibodies, we used the TSA Enhancement kit (Perkin Elmer, NEL702).

We used a Leica MZ FL III microscope (filter excitation 470/40 and emission 590LP) and a Leica DFC320 camera to view and photograph fluorescent hearts.

Embryos were fixed in 4% formaldehyde. Whole mount in situ hybridization was performed as described previously (Moorman et al., 2001).

We thank Ingeborg Hooijkaas, Rajiv Mohan and Bouke de Boer for their contributions.