Transcriptome profile of the zebrafish SAR
To localize the SAR and to distinguish it from the rest of the cardiomyocytes, the sqet33mi59BEt line was crossed with the transgenic line Tg(myl7:dsRed). In the double transgenic, dsRed-labeled CMs partially overlap with the small number of GFP-positive cells around the inflow region of the heart (Fig. 1A). To visualize the structure of the zebrafish SAR in more detail, we generated a 3D-reconstruction of the region expressing GFP in the sqet33mi59BEt transgenic line at 72 hpf (Fig. 1B) using light sheet fluorescence microscopy (LSFM). At this stage, GFP-expressing cells make a compact ring of approximately 40 cells at the inflow region of the atrium, near the sinus venosus. This location corresponded to the reported location of the zebrafish SA pacemaker [23, 26, 28]. The SAR appears to contain six-seven layers of cells at the ventral side, while only two layers of cells were observed on the dorsal side.
To obtain GFP-expressing SAR cells, intact beating hearts were manually picked at 72 hpf under a fluorescent microscope, their cells subsequently dissociated and sorted using fluorescent activated cell sorting (FACS) (Fig. 1C). The average fraction of FACS-yielded GFP-positive events obtained represented between 5% and 13% of total singlet events (Figure S1A). To assure the cells were sorted correctly, the GFP transcript levels in the two sorted fractions were compared using qPCR. GFP expression was enriched in GFP + as compared to the GFP − fraction (Figure S1B). In contrast, mRNA levels of neurogenin1 (ngn1), a neuronal gene, was not detected in GFP + and GFP − cells. Therefore, we concluded that the GFP + fraction represents the SAR cell population, while the GFP- fraction represents rest of heart (ROH).
Transcriptomic data of 19,854 genes was analyzed based on Transcript Per Million (TPM) using the oposSOM R-package . The expression portraits of SAR and ROH cell populations show two spots in opposite corners of the map, respectively, referring to opposing expression modules (Fig. 1D): in SAR, metagenes located in top-right corner of the map are overexpressed, while those in the bottom-left corner are underexpressed, and vice versa in GFP- samples, revealing antagonistic expression patterns between the two sets of samples. Pairwise correlation maps (PCM) illustrating Pearson correlation coefficients were then calculated for all replicates of the SAR and ROH. All replicates of the SAR and ROH exhibit a strong negative correlation with respect to each other, underlining the antagonistic expression patterns as seen in the sample portraits (Fig. 1E).
2. Genes and functional categories enriched in SAR cells
A total of 36 metagenes containing 730 genes were overexpressed in all 3 SAR replicates and underexpressed in ROH (Fig. 1D, top right red cluster, Table S1). In contrary, 24 metagenes and 562 genes were overexpressed in ROH and underexpressed in SAR (Fig. 1D, bottom-left red cluster, Table S1). To investigate molecular signatures specific to the SAR, we performed functional annotation analysis using oposSOM’s built-in gene-set library. Enriched GO terms among overexpressed genes in SAR cells include “cardiac muscle proliferation”, “heart contraction”, and “potassium ion transmembrane” (Fig. 2A, Table S2), which reflects the CM identity of these cells. In particular, the gene set under the “heart contraction” and “potassium ion transmembrane proteins” terms (Fig. 2B, Table S2) were overrepresented in SAR as compared to ROH. The gene encoding wnt1, a ligand of Wnt pathway, and bmp4, a ligand of BMP pathway, which have been reported to be expressed in the sinus venosus region , were among the top genes differentially expressed in SAR. Moreover, genes encoding transcription factors known to be responsible for pacemaker specification including tbx5a and tbx18 were also overexpressed (Table S1) . Unexpectedly, the pacemaker signature isl1 was not among the differently expressed genes in SAR compared to ROH. Studies have shown that the Wnt signalling promotes the specification of pacemaker cells in mammals and teleost [30, 32]. Our data set shows that genes related to the Wnt signalling were enriched in SAR as compared to ROH (Fig. 3A, top panel). We found that 55 genes of the Wnt signalling pathway were among the genes upregulated in SAR, such as frizzled homolog 9b (fzd9b), wnt2ba, wnt2bb, and wnt11. Another signalling pathway crucial for pacemaker activity is the BMP signalling [3, 33]. Our data confirms that most of the components of this pathway, including bmp2, bmp4 and bmp7b were overexpressed in SAR (Fig. 3A, lower panel).
SAN function is thought to be dependent on a 2-clock mechanism: the calcium clock and a voltage clock. The automaticity of the SAN is jointly regulated by the rhythmic spontaneous Ca2 + release from the sarcoplasmic reticulum and the cyclic activation and deactivation of membrane ion channels . The hyperpolarization-activated, cyclic nucleotide-gated channel 4 (hcn4) subunit, which is known for its role in generating the pacemaker current , was one of the SAR’s 730 overexpressed genes. In addition, calcium channel subunit, T-type channels and potassium channel related genes (kcnq1, kcnh2) known to be overexpressed in pacemaker cells , were overexpressed in the SAR transcriptome (Table S3).
In order to identify novel pacemaker-specific genes, we performed differential expression analysis, showing the two main cell type-specific groups of transcripts. At the FDR (false discovery rate) significance level of < 0.1, we identified 92 genes (49 SAR-highly-enriched genes and 43 ROH-highly-enriched) (Table S4). Interestingly, expression of some genes previously not reported in pacemaker, such as pard6a, prom2, wt1b, brd2a, atoh8, fn1b, and ldlrap1a, were higher in SAR compared to ROH (Fig. 4A, Table S4). The gene maps show the distribution of these differentially expressed genes within the original SOM, one cluster forms at the top-right corner (overexpressed in SAR) and another at the bottom-left (overexpressed in ROH) (Fig. 4B) in agreement with the sample portraits shown in Fig. 1D. The difference between these two cell populations was also evident from the Pairwise Pearson correlation heatmap (Fig. 4C).
3. Known mammalian SAN signatures are enriched in SAR
Several landmark studies have reported the transcriptome profile of mammalian SAN cells isolated by various methods. Vedantham et al. (2015) profiled the transcriptome of SAN cells isolated by laser capture microdissection from a sinus node reporter mouse line and identified differentially expressed genes between SAN and right atrium (RA) at different time points, yielding a set of SAN core genes. We found that 39 out of 96 mammalian SAN core genes identified in their study were among the top SAR-overexpressed genes. These include tbx18, hcn4, bmp4, tgfb2, and slc6a17 (Fig. 3B; Table S5). More recently, mouse and human SAN molecular signatures were reported from cells isolated from Tbx3-driven reporter mouse and human fetal SAN cells (van Eif et al., 2019). We identified 26 out of 89 conserved SAN signature genes overexpressed in zebrafish SAR (Table S5). These include cacnb1, camk4, bmp4, and hcn4. Interestingly, when we compared our SAR-specific transcriptome with results of these two studies (39 and 26 genes), we found nine common genes: shisa4, rgs6, bmp4, uchl1, sfrp5, hcn4, tub, boc, and chrnb4. These nine genes might have a crucial role in SAN development and function (Table S5). Furthermore, several known calcium ion channels (cacna1ab, slc8a4b) , extracellular matrix genes (mmp11a, thbs2b) and neuronal crest genes (foxd3,gria1a, smarca4a, and sox10), which were reported to be expressed in SAN  were also overexpressed in our SAR dataset (Table S5).
In order to assess the relevance of the zebrafish SAR transcriptome to human diseases, we ask whether any of the human orthologs of SAR-overexpressed genes were associated with congenital heart diseases. We found 24 genes implicated in 104 different ClinVar phenotype entries related to various forms of congenital heart diseases (Fig. 5, Table S7,), among which are forms of cardiac arrhythmia, including Brugada syndrome, sinus bradycardia, sudden cardiac death, and Long QT syndrome.
4. Knockdown of three SAR-overexpressed genes causes defects of heart morphology and physiology
To validate the relevance of SAR-overexpressed genes identified in our study to cardiac function, we selected three genes for functional analyses by MO knockdown. The candidate genes pard6a, prom2, and atp1a1a.2 were not previously associated with cardiac function. pard6a encodes the adapter protein involved in asymmetrical cell division and cell polarization processes in the neural tube . Prom2 localizes to basal epithelial cells and may be involved in the organization of plasma membrane microdomains . It interacts with amigo1, the known interactor of the delayed rectifier voltage-dependent potassium channel KCNB1 [40, 41]. The Atp1a1a.2-related Atp2a1a.1 was associated with cardiocytes differentiation and cell migration during heart formation . Injections of 2 ng of gene-specific, translation blocking MO resulted in developmental delay without morphological abnormalities in about 85% of larvae for pard6a, 80% for prom2, and 86% for atp1a1a.2 (Fig. 6A). Injection of 4 ng of MO resulted in a high percentage of embryos with abnormal morphology: 92.7% for pard6a, 90% for prom2, and 90.9% for atp1a1a.2 (Fig. 6B). Bright-field microscopy imaging revealed abnormalities in heart morphology. At 72 hpf, the heart in control larvae show normal cardiac looping with properly developed chambers and no cardiac edema (Fig. 6C, G). In individuals injected with 2 ng of pard6a MO, the heart morphology appeared unaffected apart from the presence of mild pericardial edema (Fig. 6D). In those injected with 2 ng of prom2 MO, both chambers were present, but the heart failed to loop and developed pericardial edema (Fig. 6E). Knockdown of atp1a1a.2 injected with 2 ng of MO did not appear to affect heart morphology compared to the control group (Fig. 6F). Injection of 4 ng of pard6a and prom2 MO resulted in a more pronounced pericardial edema, together with linearized heart tube (Fig. 6H, I). On the other hand, 4 ng of atp1a1a.2 MO had a comparatively milder phenotype: both heart chambers were distinguishable although heart looping was not complete, accompanied by mild pericardial edema (Fig. 6J). This could be due to compensation of the related atp1a1a.1.
To test if MO-induced loss-of-function affects heart physiology, heartbeat rate in MO-injected embryos were measured and compared to that of uninjected control. The average heartbeat rate for uninjected larvae at 72 hpf, was 128 beats per minute (bpm) (Fig. 6K, n = 61), which was within the reported range (Poon and Brand, 2013). Average heart rate of embryos injected with 2 ng of pard6a MO was 102.25 bpm (n = 20; p = 5.77E-04), while those injected with 4 ng of MO had a heartbeat rate of only 57 bpm (n = 20; p = 1.59E-13), which indicated severe bradycardia (Fig. 6K). Injection of prom2 MO did not result in significant changes in average heart rate (n = 20). In embryos injected with 2 ng and 4 ng of atp1a1a.2 MO, average heart rate was 111 bpm (n = 21; p = 1.26E-10) and 119 bpm (n = 20; p = 1.28E-02), suggesting bradycardia (Fig. 6K). Taken together, MO knockdown of each of the three candidate genes resulted in various forms of disruption to the heart morphology and heartbeat rate. All observed phenotypes were retained when p53 MO were co-injected with the gene-targeting MOs (Fig. 6K, S3), confirming the lack of off-target effects due to p53-dependent cell death .