Transcriptome Profile of the Sinoatrial Ring Reveals Conserved and Novel Genetic Programs of the Zebrafish Pacemaker

doi:10.21203/rs.3.rs-130687/v1

Download PDF

Research article

Transcriptome Profile of the Sinoatrial Ring Reveals Conserved and Novel Genetic Programs of the Zebrafish Pacemaker

https://doi.org/10.21203/rs.3.rs-130687/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background: Sinoatrial Node (SAN) is part of the cardiac conduction system, which controls the rhythmic contraction of the vertebrate heart. The SAN consists of a specialized pacemaker cell population that has the potential to generate electrical impulses. Although the SAN pacemaker has been extensively studied in mammalian and teleost models, including the zebrafish, their molecular nature remains inadequately comprehended.

Results: To characterize the molecular profile of the SAR and elucidate the mechanism of pacemaker function, we utilized the zebrafish transgenic line sqet33mi59BEt to isolate cells of the sinoatrial ring (SAR) of developing zebrafish embryos and profiled their transcriptome. Our analyses identified novel candidate genes and well-known conserved signaling pathways involved in pacemaker development. We show that, compared to the rest of the heart, the zebrafish SAR overexpresses several mammalian SAN pacemaker signature genes, which include hcn4 as well as those encoding calcium- and potassium-gated channels. Moreover, genes encoding components of the BMP and Wnt signaling pathways, as well as members of the Tbx family, which have previously been implicated in pacemaker development, were also overexpressed in the SAR. Among SAR-overexpressed genes, 24 had human homologues implicated in 104 different ClinVar phenotype entries related to various forms of congenital heart diseases, which suggest the relevance of our transcriptomics resource to studying human heart conditions. Finally, functional analyses of three SAR-overexpressed genes, pard6a, prom2, and atp1a1a.2, uncovered their novel role in heart development and physiology.

Conclusion: Our results established conserved aspects between zebrafish and mammalian pacemaker function and revealed novel factors implicated in maintaining cardiac rhythm. The transcriptome data generated in this study represents a unique and valuable resource for the study of pacemaker function and associated heart diseases.

Epigenetics & Genomics

zebrafish

cardiac conduction system

pacemaker

sinoatrial node

sinoatrial ring

RNA-seq

The cardiac conduction system (CCS) is an essential component of the heart across vertebrates [1, 2]. It is responsible for initiating and coordinating the electrical signals that cause rhythmic and synchronized contractions of the atria and ventricles [3]. In higher vertebrates, this system comprises of the sinoatrial node (SAN) (the primary pacemaker site), atrioventricular node (AVN) (the secondary pacemaker site), and the “wiring” of the ventricles, namely the Purkinje fibres [1]. Electrical impulses are first generated in the SAN and are rapidly propagated through the heart muscle cells (cardiomyocytes, CMs) of the atrium, resulting in atrial contraction [4]. The electrical impulses then reach the slow conducting tissues of the secondary pacemaker at the AVN, causing a delay before being passed on to the fast conducting Purkinje fibres, which form a network throughout the ventricular myocardium [3, 5]. The primary pacemaker consists of a specialized group of CMs with distinctive morphological and electrophysiological properties [6].

Studies from mammalian and fish model systems have shown that distinct components of the CCS derive from the CMs, where they were set apart early from working CMs and driven to adopt distinct morphological and physiological characteristics to perform their specialized functions [7]. The SAN, which serves as the primary pacemaker site, develops within the sinus venosus, an area located at the site of blood entry at the meeting point between the right atrium and superior vena cava. In murine embryonic heart, Tbx3-expressing CMs in the early heart tube give rise to the SAN pacemaker cells, where Tbx3 activates the pacemaker genetic program and suppresses that of working CMs [7–9]. These cells also retain the expression of LIM-homeodomain transcription factor Isl1, which is lost in differentiated CMs [10–12]. The pacemaker gene program causes the SAN pacemaker cells to retain a low proliferation rate and develop slow conduction through the expression of low conductance gap junction proteins such as Cx30 [13, 14].

The CCS is evolutionarily conserved in the building plan of the heart, and it indicates that the cellular and molecular mechanisms that drive the formation of pacemaker tissues are almost similar among vertebrates [2]. The zebrafish two-chambered heart with a single atrium and a ventricle has remarkable similarities to the human heart in terms of basal heart rate, electrophysiological properties, and action potential shape and duration [15–17]. Several zebrafish mutants have been described which display phenotypes closely resembling that found in various forms of human arrhythmia [18–22]. It suggests the high conservation of molecular pathways regulating heart conduction. However, despite its promise as a model for human cardiac arrhythmia, the zebrafish has a poorly characterized CCS. Unlike that in mammals, zebrafish pacemaker cells are difficult to distinguish morphologically from CMs. Although the zebrafish homologs of the pacemaker genes tbx2b and tbx3a are expressed in tissues which correspond to the pacemaker in mammalian heart [2, 23], no detailed analyses have been made to confirm the identity of these expression domains. Identification of the pacemaker domain in the zebrafish has relied on the mutually exclusive expression of marker gene pairs to distinguish pacemaker cells from CMs, such as isl1, which revealed a ring structure in the myocardium at the junction of sinus venosus and atrium [23]. Visualization of calcium activation wave using the fluorescently-labelled calcium transgenic line has traced pacemaker activity to a small number of cells residing within the sinoatrial ring (SAR) which served as initiation site [24]. However, to date, no live markers of the pacemaker have been available, which would enable direct visualization of its morphology and structure or isolation for molecular characterization.

A screen of enhancer trap lines generated through transposon-mediated transgenesis [25] has generated a collection of transgenic lines expressing reporter genes in different tissues or subdomains of the heart [26, 27]. In the transgenic line sqet33mi59BEt, GFP is expressed in a ring structure at the venous pole of the heart [26], which is in agreement with the known position of the pacemaker region described previously [23]. To characterize the molecular building blocks of pacemaker, we isolated cells of the SAR and profiled their transcriptome. We identified well-known conserved signaling pathways, including Wnt and BMP, as well as novel candidate genes involved in pacemaker development. Our results show that the SAR transcriptome profile exhibits similarities to that of the mammalian SAN in terms of the expression of several key genes regulating pacemaking activity as well as its development.

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 [29]. 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 [30], 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) [31]. 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 [34]. The hyperpolarization-activated, cyclic nucleotide-gated channel 4 (hcn4) subunit, which is known for its role in generating the pacemaker current [35], 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 [36], 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) [37], extracellular matrix genes (mmp11a, thbs2b) and neuronal crest genes (foxd3,gria1a, smarca4a, and sox10), which were reported to be expressed in SAN [33] 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 [38]. Prom2 localizes to basal epithelial cells and may be involved in the organization of plasma membrane microdomains [39]. 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 [42]. 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 [43].

The understanding of molecular mechanism underlying CCS development and function is one of the key challenges, which require addressing in order to develop better diagnosis and treatment for diseases related to pacemaker dysfunction. Several studies have reported the transcriptome profile of the SAN tissue isolated using laser capture microdissection (LCMD) in mouse [12] and TOMO-seq in zebrafish [30]. More recently, analysis of specific cell population of the mouse SAN was enabled by FACS to isolate specific populations of fluorescently-labelled cells from Tbx3+/Venus knock-in mice [33]. Their study revealed the SAN-enriched gene expression program, among which TBX3, SHOX2, ISL1, HOX family, BMP and NOTCH signalling, were conserved between human and mouse.

In this study we utilized a transgenic zebrafish line expressing GFP at the SAR region located between the sinus venosus and the atrium [26, 27], from which SAR cells were isolated using FACS (Figure S2). The location of GFP-positive cells in embryonic zebrafish heart and the overexpression of molecular markers shox2 and hcn4 in our data support the idea that pacemaker cells are a part of the GFP-expressing SAR and comprise cells involved in the origin of electrical current and contraction of the zebrafish heart [26] (Figure S1). Unlike the mammalian pacemaker signature, we observed that CM markers such as nppa was overexpressed while the gene encoding the slow conducting connexin Cx32.3 was underexpressed in the SAR transcriptome (Table S1). This suggests that the GFP-expressing region of the sqet33mi59BEt likely consists of a heterogenous population of cells, comprising of several cells with a definite pacemaker activity plus other cells, some of which may represent the pacemaker progenitors. This is in line with the previous observation by Chi and colleagues (2008) where the origin of activation could be traced to a small number of cells in this area.

We performed transcriptome data portraying using a self-organizing map that includes the expression patterns of all genes expressed in the SAR of 72 hpf zebrafish. Analysis of the SAR (GFP-positive cells) and ROH (GFP-negative cells) transcriptome highlighted two sets of genes: (a) genes overexpressed in SAR and (b) gene cluster overexpressed in ROH (Fig. 1). These genes are inversely expressed in both populations of cells. The 730 genes overexpressed in the SAR encompass many uncharacterized ones, indicating the potential significance of this dataset for revealing new genes and pathways involved in zebrafish CCS development and function. Interestingly, most of the SAR-highly-enriched genes have not been previously implicated in pacemaker development or function. Comparison of our dataset with recently published mouse SAN signatures [12, 33] highlights that many of the important SAN pacemaker genes are overexpressed in zebrafish SAR. However, a few known CCS markers including isl1 and shox2 were unexpectedly underrepresented in SAR as compared to ROH. This might be due to the heterogeneity of the SAR region, in which very few cells possess definite pacemaker activity and the rest maintain an intermediate character between pacemaker and working myocardium. Alternatively, other structures in the rest of the heart may also express some of these genes. Further analyses at the single cell level (currently in progress) will help us to more accurately distinguish between these different cell types within the SAR.

Loss-of-function of three novel SAR-overexpressed genes - pard6a, prom2, and atp1a1a.2, affected both heart morphology and physiology in a dose-dependent manner. Although these candidate genes have not been previously implicated in cardiac development or function, evidences exist which support their specific involvement in this process. Notably, a mutant of aPKCλ, an interactor of Pard6a, is characterized by defects in heart tube assembly resulting from a disruption to adherens junction formation [44]. Therefore, it is reasonable to believe that a similar mechanism underlies the disrupted cardiac morphology observed in pard6a knockdown. Since Pard6a morphants develop bradycardia, it seems that Pard6a plays a role in the pacemaker function, but what this exact role is remains to be established. Interestingly, bradycardia was also observed as a result of atp1a1a.2 knockdown despite its comparatively mild morphological phenotype. It is known that Na+/K+-ATPases are the major class of ion pumps in higher eukaryotes which is essential for maintaining proper heart rate by regulating intracellular ionic homeostasis of the heart cells [19, 45]. Interestingly, the ion gradient generated by Na+/K + ATPase is also known to be essential for maintaining myocardial cell junction during heart morphogenesis [46], which may partly explain the cardiac phenotype observed in atp1a1a.2 knockdown.

Thus, in this study, we characterized the molecular profile of the zebrafish SAR which exhibit pacemaker activity. The anatomical position of the SAR GFP-positive cells in transgenic line sqet33mi59BEt at the inflow region coincides with the previously reported SA pacemaker region in zebrafish heart. We show that the SAR overexpresses signature genes of mammalian SAN pacemaker, which include the hcn4 ion channel which generates the pacemaker current, as well as other calcium and potassium gated channels known to be specific to pacemaker cells. Moreover, genes encoding components of the BMP and WNT signaling pathways, as well as members of the Tbx family which have previously been implicated in pacemaker development were also overexpressed in the SAR, which suggests conservation of molecular mechanism of pacemaker development and function in zebrafish and mammals. Finally, functional analyses of three SAR-overexpressed candidate genes revealed their role in heart development and physiology, possibly through their function in the SAR. Taken together, our data represents a valuable resource for the study of pacemaker function and heart diseases associated with the defect of this structure.

Zebrafish breeding

Wild-type and transgenic zebrafish sqet33mi59BEt were maintained in the zebrafish facility of the International Institute of Molecular and Cell Biology in Warsaw (License no.PL14656251 in the register of the District Veterinary Inspectorate in Warsaw), according to standard procedures and ethical practice recommended. Embryos were grown in embryo medium at 28 °C, staged according to standard morphological criteria [47]. All experiments were done according to standard protocols established by Polish Laboratory Animal Science Association.

Imaging using Light-Sheet fluorescence and confocal microscopes

Embryos were anaesthetized at desired developmental stages and embedded into a ∼1 mm inner diameter glass capillary filled with 1.5% low-melting agarose (LMA) in embryo medium (0.03% Instant Ocean salt into double distilled water). After complete polymerization of agarose, the capillary was inserted into the sample holder, placed in a microscope chamber filled with embryo medium with 0.02% ethyl 3-aminobenzoate methanesulfonate (tricaine), and the embryo was pushed out of the capillary for imaging. The temperature of the sample chamber was maintained at 28 °C during imaging. Image acquisition was done by ZEISS Light-sheet Z.1 with W Plan-Apochromat 20x/1.0. Z-stack (thickness 3.83 µm and exposure time 60 ms). Images were saved in the LSM format and then processed using ZEN software (Zeiss). For each z-stack, maximum intensity projections were generated. For confocal microscopy, embryos were mounted in 1.5% LMA in glass-bottom dishes and after complete polymerization images were taken on the LSM800.

Embryonic Heart isolation

The zebrafish embryonic heart isolation protocol was adopted from [48]. Hearts were isolated at 72 hpf from sqet33mi59BEt zebrafish embryos. Three biological replicates of approximately 500 pooled embryos were prepared in a 5 ml low-binding tube and anaesthetized with tricaine (0.16 mg/ml in E3 medium). Embryos were washed 5x with ice-cold E3 medium and twice with ice-cold embryo destruction medium (EDM medium, L-15 medium from Gibco containing 10% FBS). Fresh EDM was added and a 5 ml syringe with 19G needle was gently pumped up and down 25 times bringing media and embryos into and out of the syringe. Media containing fragmented embryos as well as media obtained from rinsing syringe were filtered through a 100 µm nylon mesh (VWR), and flow-through was collected in a 30-mm glass petri dish. This flow-through was subsequently passed through a 40 µm nylon mesh to separate hearts from the debris and different fragments of bodies. The 40 µm mesh was inverted, and the retained material was washed off using EDM into a small petri dish. Intact and GFP positive hearts were collected using pipette under Olympus SZX16 fluorescence stereomicroscope to the low-binding tube with cold EDM medium and put on ice. The hearts were further processed using the following steps.

Heart cell dissociation and fluorescence-activated cell sorting (FACS)

GFP labelled intact hearts were chemically dissociated into a single-cell suspension for FACS. Collected hearts were centrifuged at 1500 g for 10 min at 4 °C. The pellet was dissociated with 1 ml of the solution containing 0.05% trypsin, 0.2 mM EDTA and PBS, and further incubated in this solution for 6 minutes at room temperature. 1 ml of HBSS solution containing 10 mM Hepes and 2.5% (w/v) BSA was then added to stop the reaction, and the dissociated hearts were centrifuged at 1500 g for 10 min at 4 °C. Pellet was then resuspended in 150 µl FACS Max (Genlantis). Prepared and dissociated cell suspension was passed through 0.35 µm strainer. The dissociated cells were put on the ice and sorted using BD FACSAria II. Cells were sorted into GFP positive (GFP+) and GFP negative (GFP-) fractions directly into LS TRIZOL® (Invitrogen), three biological replicates each.

RNA isolation from FACS sorted cells and assessment of RNA quantity and integrity

Total RNA was extracted from GFP + and GFP − cells using TRIzol LS® (Invitrogen) according to the manufacturer’s protocol and followed by purification using RNA Clean & Concentrator (ZYMO Research). RNA yield and quality were assessed using Agilent 2200 Tapestation with the High Sensitivity RNA Screen-Tape system according to the manufacturer’s protocol (Agilent Technologies). RNA Integrity Number (RIN^e) was in the range of between 5.8 to 6.5 for all the samples used for RNA-sEq. To obtain a high concentration, RNA was eluted in 7 µl of nuclease-free water.

cDNA synthesis for library preparation and Sequencing

cDNA synthesis was performed using SMART-Seq® v4 Ultra® Low Input RNA Kit for Sequencing (TakaRa) was used for cDNA synthesis. RNA-seq libraries were generated by Illumina Nextera XT kit according to manufacturer's protocol. Sequencing libraries were purified with Agencourt AMPure XP PCR purification beads (Beckman Coulter, USA) and fragment distribution was assessed by using High Sensitivity D1000 ScreenTape and Agilent 2200 TapeStation system (Agilent Technologies, USA). KAPA library quantification kit (Kapa Biosystems, USA) was used to quantify the libraries. RNA-seq was performed with paired-end sequencing (2 × 75 bp reads) on the NextSeq 500 (Illumina, USA). The sequencing coverage was at least 10 million reads and 4 million reads for GFP + and GFP-, respectively.

Quantitative PCR

Roche SuperScript IV First-Strand Synthesis System (Thermo Fisher Scientific) was used to obtain cDNA. Maximum amount of RNA (7 µl) was reverse transcribed to cDNA using random hexamers primers according to manufacturer’s instructions. The reaction was performed at 23ºC for 10 min, 55ºC for 10 min and finally at 80ºC for 10 min. Relative mRNA expression was quantified by using FastStart SYBR Green master mix on the LightCycler 96 instrument (Roche Life Science) with specific primer sets.

Data analysis

RNA-seq reads were mapped to the zebrafish genome assembly GRCz11 (Ensembl Genes 92). Gene expression data (TPM values; transcripts per million) were generated using Salmon [49] and further pre-processed using quantile normalization and centralization as described in [50]. Data were then used to train a self-organizing map (SOM), translating the expression profiles of all 19,854 genes measured into a two-dimensional map of 30 × 30 metagene profiles. Hereby, ‘profile’ denotes the vector of all samples’ expression values per gene or metagene. The SOM machine learning distributes the genes over the metagenes such that genes with similar profiles cluster together in the same or in closely located metagenes. After training, each metagene represents a cluster of single genes with very similar expression profiles, and the corresponding metagene profile can be interpreted as the mean profile averaged over all associated single gene profiles. The metagene expression values of each sample can be visualized by arranging them into a two-dimensional mosaic map representing the 30 × 30 metagenes, and by applying a color code presenting overexpressed metagenes in orange and red colors, intermediate ones in yellow and green, and underexpressed metagenes in light to dark blue. Each of those mosaic portraits exhibits characteristic spatial color patterns serving as a fingerprint of the transcriptional activity of the respective sample. This ‘data portraying’ approach was successfully applied to a multitude of molecular-biological data sets in the context of cancer and health research [51–53]. It is available as Bioconductor R-package ‘oposSOM’ [29]. We validated the oposSOM differential expression analysis using results of DESeq2 as an independent method analysis to compare between GFP + and ROH fractions. Using a threshold of p_adj < 0.05 and − 1.5 > = log₂FC > = 1.5, we found 858 genes upregulated and 1,059 genes down-regulated in GFP + compared to ROH (Table S6, see DESeq2_padj < 0.05). The results of DESeq2 and oposSOM largely agree with each other, with 63% (31 out of 49) of oposSOM SAR-enriched genes also identified as significantly upregulated in DESeq2 analysis.

Intersection with GWAS dataset

To intersect SAR-overexpressed genes with GWAS dataset, SNP-phenotype Associations were downloaded from ClinVar (ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/archive_2.0/2020/). The database contains 1,319,815 Variant-Trait Associations. Zebrafish orthologs were identified for all known human SNP associated genes contained in ClinVar dataset using R BioMaRt Package. A total of 37,241 unique zebrafish orthologs were found, out of which 14,008 unique SNPs corresponded to SAR-overexpressed genes and 14,524 to SAR-underexpressed genes. Further, ClinVar data were intersected with CHD specific traits, resulting in 29 unique CHD specific SNPs from SAR-overexpressed and 7 from SAR-underexpressed gene lists.

Morpholino analysis

Morpholino antisense oligomers (MO, Gene Tools, USA) were designed to block the translation of Pard6a, Prom2, and Atp1a1a.2. Two different concentrations of MO working solutions were prepared (2 ng/nl and 4 ng/nl) by diluting MO stock in deionized nuclease-free water (Polpharma, Poland) and adding 10% Phenol Red dye (Sigma Aldrich). One-cell stage embryos of ABxTL line were injected with 2 ng or 4 ng of gene-specific MO alone or together with 4 ng of p53 MO into the yolk and grown in standard conditions. For each MO type and concentration, experiment was performed with three replicates of at least 50 embryos each. Mortality, phenotype and heartbeat were assayed at 24, 48, and 72 hpf in both groups. Embryos were mounted in 2% methylcellulose and imaged from lateral side in brightfield (SZX16, Olympus). For imaging of heart morphology, 72 hpf embryos were treated with 1x tricaine for 10 minutes to stop heartbeat and mounted individually in 2% methylcellulose. Bright field images from ventral side were taken under Imager.M2 microscope (Zeiss). Heartbeat rates were measured for 1 minute under microscope by mounting each larvae in a drop of egg water on a petri dish. Experiment was performed at a stable temperature of 28 °C. Average heartbeat rate was obtained from 20 larvae for each experiment. Statistical significance between MO and control groups was determined using Student’s T-test.

Ethics approval and consent to participate

Not Applicable. According to EU Directive 2010/63/EU on the protection of animals used for scientific purposes, the earliest life-stages of animals do not fall into the regulatory frameworks of animal experimentation. This study used only zebrafish embryos less than 120 hpf which therefore does not require consent from the Local Ethics Committee.

Consent for publication

Not Applicable.

Availability of data and materials

All sequencing data have been deposited in the GEO database under accession number GSE160398.

Competing interests

The authors declare that they have no competing interests.

Funding

Research reported in this publication was supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 665778, the fellowship registration number: 2015/19/P/NZ3/03613, supported by National Science Centre (NCN), Poland, for R.M.; the Opus Grant 2016/21/B/NZ3/00354 from the National Science Centre Poland for V.K.; the Opus Grant 2018/29/B/NZ2/01010 from the National Science Centre, Poland for S.V. Microscopy, sequencing data acquisition and functional analyses was supported by the following funder: The project no. POIR.04.04.00-00-1AF0/16-00/ is carried out within the First TEAM programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund.

Authors’ contributions

R.M., V.K., and C.W. conceived the project and designed the experiments. R.M., Y.S., T.O., K.N., A.P., A.B. and L.B. performed the experiments. K.P. planned and coordinated flow cytometry/sorting experiments. H.W., H.B., R.M., C.W. and S.V. analyzed the data. The manuscript was written by R.M., Y.S., H.W., V.K., H.B. and C.W. All authors read and approved the final manuscript.

Acknowledgments

We are grateful to zebrafish core facility of the IIMCB Warsaw for excellent fish care. We thank W. Rybski and Dr. Poon K.L. for technical help. We are also obliged to Dr. R. Ramadass and other D. Stainier lab members for useful discussions during this project.

Davies F: The Conducting System of the Vertebrate Heart. Br Heart J 1942, 4(3):66–76.
Jensen B, Boukens BJ, Postma AV, Gunst QD, van den Hoff MJ, Moorman AF, Wang T, Christoffels VM: Identifying the evolutionary building blocks of the cardiac conduction system. PLoS One 2012, 7(9):e44231.
van Weerd JH, Christoffels VM: The formation and function of the cardiac conduction system. Development 2016, 143(2):197–210.
Kennedy A, Finlay DD, Guldenring D, Bond R, Moran K, McLaughlin J: The Cardiac Conduction System: Generation and Conduction of the Cardiac Impulse. Crit Care Nurs Clin North Am 2016, 28(3):269–279.
Moorman AF, de Jong F, Denyn MM, Lamers WH: Development of the cardiac conduction system. Circ Res 1998, 82(6):629–644.
Anderson RH, Yanni J, Boyett MR, Chandler NJ, Dobrzynski H: The anatomy of the cardiac conduction system. Clin Anat 2009, 22(1):99–113.
Mohan RA, Mommersteeg MTM, Dominguez JN, Choquet C, Wakker V, de Gier-de Vries C, Boink GJJ, Boukens BJ, Miquerol L, Verkerk AO et al: Embryonic Tbx3(+) cardiomyocytes form the mature cardiac conduction system by progressive fate restriction. Development 2018, 145(17).
Hoogaars WM, Barnett P, Rodriguez M, Clout DE, Moorman AF, Goding CR, Christoffels VM: TBX3 and its splice variant TBX3 + exon 2a are functionally similar. Pigment Cell Melanoma Res 2008, 21(3):379–387.
Bakker ML, Boink GJ, Boukens BJ, Verkerk AO, van den Boogaard M, den Haan AD, Hoogaars WM, Buermans HP, de Bakker JM, Seppen J et al: T-box transcription factor TBX3 reprogrammes mature cardiac myocytes into pacemaker-like cells. Cardiovasc Res 2012, 94(3):439–449.
Sun YF, Liang XQ, Najafi N, Cass M, Lin LZ, Cai CL, Chen J, Evans SM: Islet 1 is expressed in distinct cardiovascular lineages, including pacemaker and coronary vascular cells. Dev Biol 2007, 304(1):286–296.
Liang XQ, Zhang QQ, Cattaneo P, Zhuang SW, Gong XH, Spann NJ, Jiang CZ, Cao XK, Zhao XD, Zhang XL et al: Transcription factor ISL1 is essential for pacemaker development and function. J Clin Invest 2015, 125(8):3256–3268.
Vedantham V, Galang G, Evangelista M, Deo RC, Srivastava D: RNA sequencing of mouse sinoatrial node reveals an upstream regulatory role for Islet-1 in cardiac pacemaker cells. Circ Res 2015, 116(5):797–803.
Boyett MR, Inada S, Yoo S, Li J, Liu J, Tellez J, Greener ID, Honjo H, Billeter R, Lei M et al: Connexins in the sinoatrial and atrioventricular nodes. Adv Cardiol 2006, 42:175–197.
Mezzano V, Sheikh F: Cell-cell junction remodeling in the heart: possible role in cardiac conduction system function and arrhythmias? Life Sci 2012, 90(9–10):313–321.
Milan DJ, Macrae CA: Zebrafish genetic models for arrhythmia. Prog Biophys Mol Biol 2008, 98(2–3):301–308.
Leong IUS, Skinner JR, Shelling AN, Love DR: Zebrafish as a model for long QT syndrome: the evidence and the means of manipulating zebrafish gene expression. Acta Physiol 2010, 199(3):257–276.
Verkerk AO, Remme CA: Zebrafish: a novel research tool for cardiac (patho)electrophysiology and ion channel disorders. Front Physiol 2012, 3:255.
Baker K, Warren KS, Yellen G, Fishman MC: Defective "pacemaker" current (Ih) in a zebrafish mutant with a slow heart rate. Proc Natl Acad Sci U S A 1997, 94(9):4554–4559.
Shu XD, Cheng KR, Patel N, Chen FH, Joseph E, Tsai HJ, Chen JN: Na,K-ATPase is essential for embryonic heart development in the zebrafish. Development 2003, 130(25):6165–6173.
Langenbacher AD, Dong Y, Shu X, Choi J, Nicoll DA, Goldhaber JI, Philipson KD, Chen JN: Mutation in sodium-calcium exchanger 1 (NCX1) causes cardiac fibrillation in zebrafish. Proc Natl Acad Sci U S A 2005, 102(49):17699–17704.
Arnaout R, Ferrer T, Huisken J, Spitzer K, Stainier DY, Tristani-Firouzi M, Chi NC: Zebrafish model for human long QT syndrome. Proc Natl Acad Sci U S A 2007, 104(27):11316–11321.
Huttner IG, Trivedi G, Jacoby A, Mann SA, Vandenberg JI, Fatkin D: A transgenic zebrafish model of a human cardiac sodium channel mutation exhibits bradycardia, conduction-system abnormalities and early death. J Mol Cell Cardiol 2013, 61:123–132.
Tessadori F, van Weerd JH, Burkhard SB, Verkerk AO, de Pater E, Boukens BJ, Vink A, Christoffels VM, Bakkers J: Identification and functional characterization of cardiac pacemaker cells in zebrafish. PLoS One 2012, 7(10):e47644.
Weber M, Scherf N, Meyer AM, Panakova D, Kohl P, Huisken J: Cell-accurate optical mapping across the entire developing heart. Elife 2017, 6.
Parinov S, Kondrichin I, Korzh V, Emelyanov A: Tol2 transposon-mediated enhancer trap to identify developmentally regulated zebrafish genes in vivo. Dev Dyn 2004, 231(2):449–459.
Poon KL, Liebling M, Kondrychyn I, Brand T, Korzh V: Development of the cardiac conduction system in zebrafish. Gene Expr Patterns 2016, 21(2):89–96.
Poon KL, Liebling M, Kondrychyn I, Garcia-Lecea M, Korzh V: Zebrafish cardiac enhancer trap lines: new tools for in vivo studies of cardiovascular development and disease. Dev Dyn 2010, 239(3):914–926.
Chi NC, Shaw RM, Jungblut B, Huisken J, Ferrer T, Arnaout R, Scott I, Beis D, Xiao T, Baier H et al: Genetic and physiologic dissection of the vertebrate cardiac conduction system. PLoS Biol 2008, 6(5):e109.
Loffler-Wirth H, Kalcher M, Binder H: oposSOM: R-package for high-dimensional portraying of genome-wide expression landscapes on bioconductor. Bioinformatics 2015, 31(19):3225–3227.
Burkhard SB, Bakkers J: Spatially resolved RNA-sequencing of the embryonic heart identifies a role for Wnt/beta-catenin signaling in autonomic control of heart rate. Elife 2018, 7.
Burczyk MS, Burkhalter MD, Tena TC, Grisanti LA, Kauk M, Matysik S, Donow C, Kustermann M, Rothe M, Cui Y et al: Muscarinic receptors promote pacemaker fate at the expense of secondary conduction system tissue in zebrafish. JCI Insight 2019, 4(20).
Bressan M, Liu G, Mikawa T: Early mesodermal cues assign avian cardiac pacemaker fate potential in a tertiary heart field. Science 2013, 340(6133):744–748.
van Eif VWW, Stefanovic S, van Duijvenboden K, Bakker M, Wakker V, de Gier-de Vries C, Zaffran S, Verkerk AO, Boukens BJ, Christoffels VM: Transcriptome analysis of mouse and human sinoatrial node cells reveals a conserved genetic program. Development 2019, 146(8).
Joung B, Ogawa M, Lin SF, Chen PS: The calcium and voltage clocks in sinoatrial node automaticity. Korean Circ J 2009, 39(6):217–222.
Stieber J, Herrmann S, Feil S, Loster J, Feil R, Biel M, Hofmann F, Ludwig A: The hyperpolarization-activated channel HCN4 is required for the generation of pacemaker action potentials in the embryonic heart. Proc Natl Acad Sci U S A 2003, 100(25):15235–15240.
Grant AO: Cardiac ion channels. Circ Arrhythm Electrophysiol 2009, 2(2):185–194.
Tyagi S, Ribera AB, Bannister RA: Zebrafish as a Model System for the Study of Severe CaV2.1 (alpha1A) Channelopathies. Front Mol Neurosci 2019, 12:329.
Chen X, An Y, Gao Y, Guo L, Rui L, Xie H, Sun M, Lam Hung S, Sheng X, Zou J et al: Rare Deleterious PARD3 Variants in the aPKC-Binding Region are Implicated in the Pathogenesis of Human Cranial Neural Tube Defects Via Disrupting Apical Tight Junction Formation. Hum Mutat 2017, 38(4):378–389.
Fargeas CA, Florek M, Huttner WB, Corbeil D: Characterization of prominin-2, a new member of the prominin family of pentaspan membrane glycoproteins. J Biol Chem 2003, 278(10):8586–8596.
Peltola MA, Kuja-Panula J, Liuhanen J, Voikar V, Piepponen P, Hiekkalinna T, Taira T, Lauri SE, Suvisaari J, Kulesskaya N et al: AMIGO-Kv2.1 Potassium Channel Complex Is Associated With Schizophrenia-Related Phenotypes. Schizophr Bull 2016, 42(1):191–201.
Luck K, Kim DK, Lambourne L, Spirohn K, Begg BE, Bian W, Brignall R, Cafarelli T, Campos-Laborie FJ, Charloteaux B et al: A reference map of the human binary protein interactome. Nature 2020, 580(7803):402–408.
Pott A, Bock S, Berger IM, Frese K, Dahme T, Kessler M, Rinne S, Decher N, Just S, Rottbauer W: Mutation of the Na(+)/K(+)-ATPase Atp1a1a.1 causes QT interval prolongation and bradycardia in zebrafish. J Mol Cell Cardiol 2018, 120:42–52.
Robu ME, Larson JD, Nasevicius A, Beiraghi S, Brenner C, Farber SA, Ekker SC: p53 activation by knockdown technologies. PLoS Genet 2007, 3(5):e78.
Horne-Badovinac S, Lin D, Waldron S, Schwarz M, Mbamalu G, Pawson T, Jan Y, Stainier DY, Abdelilah-Seyfried S: Positional cloning of heart and soul reveals multiple roles for PKC lambda in zebrafish organogenesis. Curr Biol 2001, 11(19):1492–1502.
Schwinger RH, Bundgaard H, Muller-Ehmsen J, Kjeldsen K: The Na, K-ATPase in the failing human heart. Cardiovasc Res 2003, 57(4):913–920.
Cibrian-Uhalte E, Langenbacher A, Shu X, Chen JN, Abdelilah-Seyfried S: Involvement of zebrafish Na+,K + ATPase in myocardial cell junction maintenance. J Cell Biol 2007, 176(2):223–230.
Kimmel CB, Ballard WW, Kimmel SR, Ullmann B, Schilling TF: Stages of embryonic development of the zebrafish. Dev Dyn 1995, 203(3):253–310.
Burns CG, MacRae CA: Purification of hearts from zebrafish embryos. Biotechniques 2006, 40(3):274, 276, 278 passim.
Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C: Salmon provides fast and bias-aware quantification of transcript expression. Nat Methods 2017, 14(4):417–419.
Wirth H, Loffler M, von Bergen M, Binder H: Expression cartography of human tissues using self organizing maps. BMC Bioinformatics 2011, 12:306.
Kunz M, Loffler-Wirth H, Dannemann M, Willscher E, Doose G, Kelso J, Kottek T, Nickel B, Hopp L, Landsberg J et al: RNA-seq analysis identifies different transcriptomic types and developmental trajectories of primary melanomas. Oncogene 2018, 37(47):6136–6151.
Binder H, Willscher E, Loeffler-Wirth H, Hopp L, Jones DTW, Pfister SM, Kreuz M, Gramatzki D, Fortenbacher E, Hentschel B et al: DNA methylation, transcriptome and genetic copy number signatures of diffuse cerebral WHO grade II/III gliomas resolve cancer heterogeneity and development. Acta Neuropathol Commun 2019, 7(1):59.
Loeffler-Wirth H, Kreuz M, Hopp L, Arakelyan A, Haake A, Cogliatti SB, Feller AC, Hansmann ML, Lenze D, Moller P et al: A modular transcriptome map of mature B cell lymphomas. Genome Med 2019, 11(1):27.

ARRIVEAuthorChecklistFull1.pdf
FigureS1.pdf
Figure S1
FACS sorting and qPCR validation of GFP positive and GFP negative samples
(A) FACS sorting of the embryonic hearts at 72hpf shows distinct GFP positive and GFP negative cell populations. (B) qPCR validation of the samples was done to check for GFP and neurogenin mRNA levels in the GFP positive and GFP negative samples. Y-axis represents fold change in percentage. qPCR was performed using ef-1α as endogenous control. RNA was extracted, reverse transcribed to cDNA and qRT-PCR reactions run using gene-specific primers to analyse mRNA levels.
FigureS2.pdf
Figure S2
Comparison with mammalian SAN signatures
Heat map showing conserved mammalian SAN genes reported in [33] that have been compared with our SAR and ROH data set. 26 out of 39 genes among the conserved mammalian signature genes were also up-regulated in our dataset.
FigureS3revised.pdf
Figure S3 Comparison of morpholino injections Comparison between those injected with only 4ng of gene-specific MO and those co-injected with p53 MO revealed no significant difference in morphological phenotype.
TableS1.xlsx
Table S1: - List of genes under expressed in SAR or overexpressed in ROH
TableS2.xlsx
Table S2: - List of differential GO sets, and their associated genes overexpressed among SAR (p < 0.1)
TableS3.xlsx
Table S3: - List of genes within the GO categories Calcium ion transport and potassium ion transport
TableS4.xlsx
Table S4: - SAR-highly-enriched genes
TableS5.xlsx
Table S5: - Comparison of SAR-overexpressed genes with mammalian SAN signatures
TableS6.xlsx
Table S6: - Intersection between OposSOM and DESeq2 overexpressed genes (padj <0.05)
TableS7.xlsx
Table S7: - List of SAN overexpressed genes associated with ClinVar phenotype

Download PDF

Editorial decision: Major revision
27 Jan, 2021
Review #1 received at journal
17 Jan, 2021
Review #2 received at journal
05 Jan, 2021
Reviewer #2 agreed at journal
29 Dec, 2020
Reviewers invited by journal
26 Dec, 2020
Reviewer #1 agreed at journal
26 Dec, 2020
Editor assigned by journal
19 Dec, 2020
Submission checks completed at journal
17 Dec, 2020
Editor invited by journal
14 Dec, 2020

You are reading this latest preprint version

Transcriptome Profile of the Sinoatrial Ring Reveals Conserved and Novel Genetic Programs of the Zebrafish Pacemaker

Status:

Version 1

Abstract

Figures

Background

Results

Transcriptome profile of the zebrafish SAR

2. Genes and functional categories enriched in SAR cells

3. Known mammalian SAN signatures are enriched in SAR

4. Knockdown of three SAR-overexpressed genes causes defects of heart morphology and physiology

Discussion

Conclusions

Methods

Zebrafish breeding

Imaging using Light-Sheet fluorescence and confocal microscopes

Embryonic Heart isolation

Heart cell dissociation and fluorescence-activated cell sorting (FACS)

RNA isolation from FACS sorted cells and assessment of RNA quantity and integrity

cDNA synthesis for library preparation and Sequencing

Quantitative PCR

Data analysis

Intersection with GWAS dataset

Morpholino analysis

Declarations

References

Supplementary Files

Status:

Version 1