Transcriptional and Epigenetic Landscape of Cardiac Pacemaker Cells: Insights Into Cellular Specialization in the Sinoatrial Node

Cardiac pacemaker cells differentiate and functionally specialize early in embryonic development through activation of critical gene regulatory networks. In general, cellular specification and differentiation require that combinations of cell type-specific transcriptional regulators activate expression of key effector genes by binding to DNA regulatory elements including enhancers and promoters. However, because genomic DNA is tightly packaged by histones that must be covalently modified in order to render DNA regulatory elements and promoters accessible for transcription, the process of development and differentiation is intimately connected to the epigenetic regulation of chromatin accessibility. Although the difficulty of obtaining sufficient quantities of pure populations of pacemaker cells has limited progress in this field, the advent of low-input genomic technologies has the potential to catalyze a rapid growth of knowledge in this important area. The goal of this review is to outline the key transcriptional networks that control pacemaker cell development, with particular attention to our emerging understanding of how chromatin accessibility is modified and regulated during pacemaker cell differentiation. In addition, we will discuss the relevance of these findings to adult sinus node function, sinus node diseases, and origins of genetic variation in heart rhythm. Lastly, we will outline the current challenges facing this field and promising directions for future investigation.

robust electrical automaticity with weak intercellular electrotonic coupling (Figure 1). In the context of the SAN, these cellular specializations allow pacemaker cells to perform their primary function -dynamic impulse generation at rates that ensure cardiac output meets metabolic need.
Over the past several decades, a prevailing model for the cellular basis of pacemaker cell function has emerged in which mutual entrainment of a membrane clock (driven by hyperpolarization-activated cation current) and a calcium clock (driven by spontaneous sub-sarcolemmal calcium release) confers robust automatic firing (Lakatta et al., 2010). In parallel, investigations at the tissue level have shown that the SAN has a distributed architecture that spans the region from the superior vena cava-right atrium junction to the inferior vena cava, bounded by the crista terminalis and the interatrial septum (Boyett et al., 2000). Nodal tissue consists of a fibrous matrix in which centrally located pacemaker cells are poorly electrotonically coupled with the atrial myocardium while pacemaker cells at the SAN periphery, where impulse transmission occurs, are more robustly coupled to the atrium (Joyner and van Capelle, 1986;De Maziere et al., 1992;Nikolaidou et al., 2012). This architecture allows for source-sink matching between the SAN and atrium without compromising automaticity.
Depending on neurohormonal state and electrolyte concentration, the nodal impulse can arise from different locations within the intercaval groove, generating an electrical wavefront that ultimately exits the SAN through species-specific pathways (Nikolaidou et al., 2012). Aligned with this conceptual framework, recent evidence has demonstrated the existence of discrete nodal structures within the SAN region that control heart rate under different conditions, with a superior SAN (located at the superior vena cava-RA junction) containing the dominant pacemaker at rapid heart rates and an inferior SAN (located near the inferior vena cava) operating at lower rates (Brennan et al., 2020). Finally, the atrioventricular node, which can function as a subsidiary pacemaker, also contains pacemaker-like cells that exhibit automaticity and hyperpolarization-activated cation current.
Remarkably, despite the anatomical and physiological complexity in heart rate regulation, both the cellular physiology of pacemaker cells and the tissue architecture in the SAN are broadly similar among diverse mammalian species despite dramatic variation in body mass and heart rate, suggesting that the determinants of cardiac automaticity constitute a robust and deeply evolutionarily conserved morphogenetic and cellular program (Bleeker et al., 1980;Opthof et al., 1985). In keeping with this notion, many of the specific molecular components required for pacemaker cell automaticity in the SAN have been identified and are present in pacemaker cells isolated from widely divergent vertebrate species (fish, mouse, and human), including hyperpolarization-activated, cyclic-nucleotide gated ion channels (Hcn4, Hcn1) and lower conductance gap junctions (connexin-45, connexin 30.2) (Blaschke et al., 2007;Tessadori et al., 2012). Conversely, high-conductance gap junctions such as connexin-43 and connexin-40, as well as channels associated with rapid conduction (NaV1.5) and hyperpolarized resting potential (Kir2.1), are present at much lower levels in pacemaker cells than in atrial cardiomyocytes (Mommersteeg et al., 2007b). More recently, genome wide expression profiles using RNA sequencing on mouse and human SAN have identified numerous common transcripts that are differentially expressed between SAN and non-SAN cardiomyocytes in both species, indicative of a shared transcriptional program that underlies the phenotypic similarities among mammalian pacemaker cells in different species (Figure 1; Vedantham et al., 2015;Goodyer et al., 2019;van Eif et al., 2019).
It is thus clear that (1) the problem of generating a robust mechanism for heart rate regulation was solved early in vertebrate evolution through the differentiation, proliferation, and architectural assembly of a specialized subtype of cardiomyocyte at the venous pole of the heart; and (2) pacemaker cells become functionally specialized by executing a conserved gene expression program with increased transcription of genes responsible for spontaneous firing and autonomic connectivity, and decreased transcription of genes associated with force generation and rapid conduction.
By connecting cellular specialization at the phenotypic level to cell type-specific gene expression programs, this unified model for the origin of vertebrate heart rhythm raises several important questions that will be the focus of the present article. First, what are the upstream regulators that allow pacemaker cells to specialize and differentiate during heart development? Second, how do these regulators turn on the SAN gene program and how is this program maintained? Third, what is the relevance of these regulatory pathways to human genetic variation and human disease? Historically, our ability to address these questions has been limited by the lack of genetic tools to mark and track pacemaker cells in model systems, as well as the difficulty of isolating sufficient quantities of purified pacemaker cells for molecular analysis. Therefore, in the present work, we will also touch on technological breakthroughs in the past decade that have allowed us to define gene expression and chromatin dynamics in pacemaker cells with unprecedented genomic resolution.

A DISTINCT SET OF TRANSCRIPTION FACTORS IS ENRICHED IN PACEMAKER CELLS AND IS REQUIRED TO ESTABLISH A SPECIALIZED GENE PROGRAM
In general, tissue-specific gene expression is accomplished through cell-type specific transcription factors (TFs) that bind combinatorially and often synergistically to cis-acting genomic DNA sequences known as enhancers that can be located far from their target genes (more on this below). These binding events, in turn, permit the formation of transcriptional complexes between TF-bound enhancers and basal promoters of target genes to activate expression. The specific complement of activated transcriptional regulators available to bind genomic DNA has been termed the "nuclear regulatory environment." The FIGURE 1 | Differences between pacemaker cells and working cardiomyocytes. (A) Characteristics of typical pacemaker cells: (top panel) micrograph of an isolated pacemaker cell taken from an adult Hcn4-GFP mouse demonstrates an elongated spindle-like morphology; (middle panel) a typical pacemaker cell action potential with elevated resting membrane potential and diastolic depolarization; (bottom panel), a list of key functional proteins enriched in pacemaker cells. (B) For comparison, typical ventricular cardiomyocyte characteristics are shown, including (top panel) a micrograph of an adult mouse ventricular cardiomyocyte stained for alpha-actinin (green) to demonstrate sarcomeres and connexin-43 (red) to demonstrate gap junctions; (middle panel) a typical working cardiomyocyte is quiescent with a negative resting potential, rapid upstroke, and plateau phase; (bottom panel) key functional genes expressed in working cardiomyocytes that are reduced in pacemaker cells.
nuclear regulatory environment in a given cell, in combination with the distribution of available enhancers, determines which genes are expressed.
Pacemaker cells express a set of transcriptional regulators that are absent in atrial cardiomyocytes: Tbx3, Tbx18, Shox2, and Isl1 (Hoogaars et al., 2004;Christoffels et al., 2006;Blaschke et al., 2007;Sun et al., 2007;Weinberger et al., 2012). These factors could therefore function as cell-intrinsic regulators that directly activate gene expression of pacemaker cell-specific genes and/or repress expression of the atrial gene program. Conversely, there are transcriptional regulators that are present in atrial cardiomyocytes that are excluded from pacemaker cells, including Nkx2.5 and Pitx2c, that may be important for activating atrial genes and repressing pacemaker cell genes (Christoffels et al., 2006;Mommersteeg et al., 2007b). Over the past 15 years, several rigorous in vivo global and conditional loss of function studies have been carried out in animal models with each of these factors, leading to a general model for transcriptional control summarized in this section.

Repression of the Atrial Gene Program in Pacemaker Cells
Several pacemaker cell-enriched transcription factors directly repress expression of atrial-enriched genes in the SAN. By blocking differentiation toward a working cardiomyocyte cell fate, this inhibition is hypothesized to permit differentiation along the path to pacemaker cells.
Tbx3, a T-box transcription factor whose expression within the SAN is conserved across several species, is the most wellcharacterized factor to play this role. Tbx3 binds to the regulatory elements that control expression of working myocardial genes such as Nppa and Gja5, possibly competing with Tbx5 (an activating T-box factor) to inhibit their expression and thereby inhibit chamber differentiation within the SAN (Hoogaars et al., 2004). Mice lacking Tbx3 demonstrate ectopic expression of atrial genes within the SAN with SAN hypoplasia, and exhibit embryonic lethality, whereas conditional deletion of Tbx3 within the adult SAN also results in ectopic atrial gene expression in the SAN and sinus arrhythmias (Hoogaars et al., 2007). Underscoring its primary role as a repressor, transgenic misexpression of Tbx3 within atrial myocardium can repress the atrial gene program and even result in some ectopic arrhythmogenic foci, but it does not directly activate expression of pacemaker cell genes ectopically.
Analogous to Tbx3, Tbx18 also appears to function as a repressor of atrial genes such as Gja5 (Wiese et al., 2009;Kapoor et al., 2011). However, as the central sinus node fails to form in the Tbx18 −/− mouse, additional specific targets of Tbx18 within the SAN have not been defined using loss of function models (Christoffels et al., 2006). In a large animal model, Tbx18 can reprogram ventricular myocardium to exhibit automaticity and transiently pace the heart, although the transcriptional mechanisms underlying this effect have not been determined (Kapoor et al., 2011(Kapoor et al., , 2013. The homeodomain transcription factor Shox2 also restricts atrial gene expression within pacemaker cells to promote SAN development. Complete loss of Shox2 causes embryonic lethality with hypoplasia of the sinus venosus, reduced embryonic heart rate, and ectopic expression of atrial genes in the SAN domain (Blaschke et al., 2007;Espinoza-Lewis et al., 2009). Recent findings have suggested a competitive mechanism whereby Shox2 directly represses Nkx2.5 expression in pacemaker cells and occupies sites normally bound by Nkx2.5 in atrial myocardium to repress atrial gene expression, thereby maintaining indirectly the expression of genes such as Hcn4, Tbx3, and Isl1 within the SAN (Ye et al., 2015). In myocardium at the interface between the SAN and the atrium, where Shox2 and Nkx2.5 are both expressed at higher levels, this competition may produce an intermediate cellular phenotype that facilitates electrophysiological sourcesink matching between the SAN interior and the atrium .

Repression of Pacemaker Cell Program in Atrial Cardiomyocytes
Within the working myocardium, the cardiac TF Nkx2.5 may act to repress expression of pacemaker cell genes. For example, conditional deletion of Nkx2.5 in ventricular myocardium results in upregulation of Hcn1, a gene normally restricted to pacemaker cells (Pashmforoush et al., 2004). Further supporting such a role, Nkx2.5 is among the few core cardiac TFs that is downregulated in pacemaker cells and their embryonic precursors in the sinus venosus (Christoffels et al., 2006). Transgenic misexpression of Nkx2.5 in pacemaker cells leads to repression of both Hcn4 and Tbx3, key components of the pacemaker cell-specific gene regulatory network (Espinoza-Lewis et al., 2011).
The homeodomain transcription factor Pitx2c is expressed throughout the left atrial and pulmonary vein myocardium, where it is largely responsible for encoding left sided identity within the atrium (Liu et al., 2002;Mommersteeg et al., 2007a). Accordingly, Pitx2c −/− mice can develop right atrial isomerism and form a secondary SAN in the left atrium that is normally repressed by the presence of Pitx2c (Mommersteeg et al., 2007b). Importantly, the Pitx2 locus harbors genomic variation that contributes to inherited susceptibility to atrial fibrillation, an arrhythmia that is triggered by focal firing from the region of the pulmonary venous myocardium and left-sided SAN remnant (Roselli et al., 2018;Zhang et al., 2019).

Activation of Specialized Gene Expression Programs in Pacemaker Cells and Atrial Myocardium
One of the major remaining unknowns in pacemaker cell transcriptional regulation is how pacemaker cell-specific gene expression is activated in the SAN. Data from avians and zebrafish have suggested that pacemaker cell progenitors arise from mesodermal cells located at the periphery of the embryonic heart fields in response to a Wnt signal (Bressan et al., 2013;Ren et al., 2019), and subsequently adopt their unique phenotype and transcriptional program as they become integrated with the atrial myocardium in a dynamic cellular process that occurs after heart looping (Bressan et al., 2018;Thomas et al., 2021).
Once committed, pacemaker cells express the broad cardiac transcription factors Gata4, Gata6, Mef2C, and Tbx5 at relatively high levels during development, similar to atrial cardiomyocytes (Vedantham et al., 2015;van Eif et al., 2019). These factors cooperatively activate the cardiac gene expression program, and their binding sites are overrepresented near both atrial cardiomyocyte and pacemaker cell genes (Galang et al., 2020;van Eif et al., 2020). Accordingly, depletion of any of these factors profoundly affects SAN development, heart rhythm, and heart development more broadly (Molkentin et al., 1997;Bruneau et al., 2001;Moskowitz et al., 2007;Kuratomi et al., 2009). Thus, although these factors are undoubtedly relevant to directing pacemaker cell-specific genes to the SAN, they do not fully account for the cell type-specificity of expression seen in pacemaker cells. One possibility is that within pacemaker cells, the repressive TFs Tbx3 and Shox2 bind to loci important for working myocardium and thereby cause Mef2C, Gata4, and Tbx5 to redistribute occupancy to pacemaker cell-specific loci where they activate transcription. This mechanism would not necessarily require any additional pacemaker cell-specific factors to explain pacemaker cell-specific gene expression.
More recently, the LIM homeodomain TF Isl1 has emerged as a potential activator of pacemaker cell-specific gene expression that could explain how pacemaker cells differentiate along a transcriptional path that is distinct from the atrial cardiomyocyte program. Although Isl1 is expressed broadly in cardiac progenitor cells of the second heart field, its expression becomes restricted to the SAN and its precursors by mid-development onward (Cai et al., 2003;Sun et al., 2007). Conditional deletion of Isl1 within the SAN results in dysregulated gene expression with major downregulation of pacemaker cell-specific genes and some upregulation of atrial genes (Liang et al., 2015). Moreover, these Isl1 SAN-conditional loss of function embryos exhibit sinus bradycardia and SAN hypoplasia, confirming a critical role for this factor in function and morphogenesis of the SAN. Analysis of promoter-adjacent regions of SAN genes downregulated after deletion of Isl1 demonstrated enrichment of Isl1 binding sites, a finding that raises the likelihood of an activating role for Isl1 in pacemaker cell-specific gene expression (Vedantham et al., 2015).

REGULATORY BASIS FOR CELLULAR SPECIALIZATION: ROLE OF GENOME ARCHITECTURE AND CHROMATIN ACCESSIBILITY
While the genetic models discussed in the previous section have provided insight into the factors that likely activate and repress cell-type specific gene expression patterns in pacemaker cells, the mechanisms that connect these factors to their target genes have not been elucidated. To a great extent, this knowledge gap has resulted from the inability to identify specific genomic loci (enhancers) where combinations of pacemaker cell and other TFs bind directly and regulate transcription of specific genes. In the absence of such direct binding information, the connection between TFs and their targets cannot be definitively established, and a common logic to combinatorial activation of pacemaker cell gene expression cannot be defined. Thus, identifying enhancers that exhibit pacemaker cell-specific activity would add another crucial piece to this complex puzzle. Before reviewing recent experiments that have tackled this problem, this section briefly outlines our current model for how enhancers and promoters interact with each other in the context of the nucleus, and the experimental methods that can be used to identify enhancers.

Organization of Genomic DNA in the Nucleus Facilitates Transcription
In mammalian cells, chromosomes each have distinct territories within the nucleus, and within each territory, transcriptionally inactive genomic loci are sequestered at the nuclear periphery, while actively transcribed regions are positioned in the interior of the nucleus (Figure 2A; van Steensel and Belmont, 2017). These transcriptionally active chromatin regions are organized into topologically associated domains (TADs), which can range in length from kilobases to megabases and can contain one or many genes along with non-coding regions that contain enhancers ( Figure 2B; Dixon et al., 2016). TADs are formed when the multifunctional cohesin protein complex is recruited by chromatin-bound CCCTC-binding factor (CTCF). Cohesin then extrudes chromatin through the center of its ring-like structure to form a loop. Loop extrusion occurs until Cohesin encounters CTCF bound in the inverted orientation, thereby fixing the boundaries of the TAD along the genome and defining the regions for potential interaction. Importantly, during the process of loop extrusion, distant chromatin elements along the TAD are brought into close proximity, permitting contact of TF-bound enhancers with target promoters, and thereby initiating transcription ( Figure 2B; Hanssen et al., 2017). Genomic regions such as enhancers and promoters within a TAD can interact whereas interactions across TAD boundaries rarely occur.
In order for DNA binding TFs to access enhancers, histone proteins undergo covalent modifications at specific sites that relax their otherwise tight association with genomic DNA, which then assumes an "open chromatin" conformation ( Figure 2C; Boland et al., 2014). These histone modifications are regulated by specific classes of enzymes, and in some cases, by a group of TFs known as "pioneer factors" that can bind to their sites in nucleosomal DNA in a "closed configuration" and recruit factors that result in opening of chromatin and binding of larger transcriptional complexes (Larson et al., 2021). In this manner, enhancers can "read out" the nuclear regulatory environment to determine whether a specific gene should be transcribed. Thus, in order to understand how transcription factors regulate their target genes in differentiated cell types, comparing the chromatin landscape among closely related cell types can uncover the specific set of enhancers that activate tissue-specific gene expression.

Experimental Approaches to Epigenetic Profiling: Relevance to Pacemaker Cells
Because chromatin at transcriptionally active enhancers must be accessible to DNA binding factors, the Assay for Transposase Accessible Chromatin with Sequencing (ATAC-seq) offers a conceptually straightforward method to identify regulatory elements genome-wide. In this technique, a bacterial Tn5 transposase is used to insert oligonucleotide tags into regions of accessible chromatin (locations where histones are in the relaxed or open conformation). Sequencing these regions in a given cell type allows for high-resolution determination of regions of genome accessibility. ATAC-seq datasets can be compared among closely related cell types to identify loci of differential accessibility that might contain tissue-specific enhancers ( Figure 2D; Buenrostro et al., 2013). ATAC-seq can be performed successfully on tens of thousands of cells and more recently, single cell ATAC-seq methodologies are available (Buenrostro et al., 2015), making this a practical tool for use with cardiac pacemaker cells directly isolated from SAN tissue.
Several factors complicate the interpretation of ATAC-seq data: First, differentially accessible genome regions are not always transcriptionally active enhancers; Second, a single enhancer can have more than one gene target; Third, because most TADs contain multiple genes, it is not always clear how to assign a putative enhancer to one or more specific gene targets within its TAD. One possible experimental approach to address these problem is to define contact frequencies among distant genomic loci through chromosome conformation capture and its variants (Sati and Cavalli, 2017). By covalently crosslinking genomic DNA and then sequencing crosslinked fragments, libraries can be generated in which sequences include fragments from an enhancer and the distant promoter it regulates since these will be in close spatial proximity during active transcription. With enough library complexity and sequencing depth, a genomewide proximity map at tens of kilobase resolution can be readily generated to define mutually interacting regions within each TAD ( Figure 2E). Variants of this approach, including promotercapture Hi-C, can enrich libraries for promoter-containing fragments, thereby compiling a high-fidelity list of promoterenhancer interactions in a given cell type Schoenfelder et al., 2015;Montefiori et al., 2018).
Unfortunately, the limiting number of pacemaker cells (∼10,000 per heart) presents an insurmountable barrier for chromosome conformation capture, which ordinarily requires 10-20 million cells. In the future, scalable in vitro strategies to differentiate pacemaker cells from induced pluripotent stem cells may provide enough material, or, alternatively, lower input chromosome conformation techniques may be become available. Until then, connecting putative enhancers in pacemaker cells with target genes will require experimental validation using in vivo models. As detailed below, these three datasets exhibited a remarkable congruence and have thereby provided novel insights into a deeply conserved cis-regulatory architecture in pacemaker cells ( Figure 3A).

Motif Enrichment Analysis and in vivo Transgenic Assays Demonstrate a Deeply Conserved Regulatory Code in Pacemaker Cells
Because regions of accessible chromatin are bound by tissue-specific transcription factors, binding sites for lineagedetermining TFs should be overrepresented in chromatin regions that are differentially accessible between closely related cell types. Accordingly, both Galang et al. and van Eif et al. looked at DNA binding motifs enriched in differentially accessible regions as compared to background sequences using an unbiased search algorithm. The results revealed a striking convergence in the motifs in both sets of differentially accessible regions. Unsurprisingly, known cardiomyocyte transcription factor motifs including Gata and T-box motifs were overrepresented. More strikingly, differentially accessible regions from both data sets, drawn from very different types of samples (mouse primary pacemaker cells versus human iPSC-derived SANLPC), showed overrepresentation of Isl1 TF binding motifs, as well as binding motifs for Meis1, a TALE-class homeodomain factor not previously studied in the context of SAN gene expression ( Figure 3B). Similar findings were observed in the data set from Fernandez-Perez et al. (2019). The finding of overrepresentation of Isl1 motifs among the differentially accessible peaks in pacemaker cells, in particular, provides strong evidence that Isl1 is a transcriptional activator in pacemaker cells and that, together with interacting partners, is likely to be a key player in activating the genes that endow pacemaker cells with their unique phenotypes.
To test whether differentially accessible ATAC-seq peaks were sufficient to function as regulatory elements in vivo, both Galang et al. and van Eif et al. cloned selected differentially accessible ATAC-seq peaks they identified and tested them using enhancer-reporter transgenesis in mice and in zebrafish. In Galang et al.,17 differentially accessible peaks were selected and tested using EnSERT, a locus-specific transgenesis assay, or using an hsp68LacZ reporter relying on random insertion. Of these, 4 ATAC-seq peaks directed reporter activity to the SAN primordium in mouse embryos; notably, these peaks were located within TADs that encompassed Hcn4, Rgs6, Ptgfr, and Isl1, all of which are differentially expressed in pacemaker cells. The enhancer at the Isl1 locus, in particular, exhibited remarkable specificity for pacemaker cells in the developing and postnatal heart, suggesting that it functions within the pacemaker cell lineage from a very early developmental stage onward.
In a similar vein, van Eif et al. tested 8 human differentially accessible ATAC-seq peaks located within TADs encompassing SHOX2, ISL1, and TBX3, and identified several with expression in the heart and sinus venosus in zebrafish, of which 3 had activity in the sinus venosus region of mice (2 at the SHOX2 locus and one at the TBX3 locus). Notably, the human region syntenic to the Isl1 enhancer identified in mouse pacemaker cells by Galang et al. was also strongly differentially accessible in the human ATAC-seq data from van Eif et al. providing a convincing example of cross-species conservation of genomic regulatory architecture. Indeed, both studies found that enhancer function was conserved across species to a remarkable extenthuman sequences were active in mouse and fish, and mouse sequences were active in fish -often with the expected tissue-specific expression pattern. The fact that the nuclear regulatory environment in zebrafish pacemaker cells can "read out" both human and mouse DNA to achieve tissue-specific expression indicates that the regulatory kernel involving Isl1, Shox2, and other pacemaker cell TFs is an ancient vertebrate regulatory module that has been conserved over hundreds of millions of years.

Enhancer Knockouts Demonstrate Required Roles for Enhancers in Gene Regulation, Development, and Function
To test the functional relevance of these enhancers, both Galang et al. and van Eif et al. deleted regulatory sequences from the FIGURE 4 | Sinoatrial node hypoplasia and dysfunction after deletion of a key pacemaker cell-specific enhancer. (A) Three-dimensional reconstruction of light sheet fluorescence image of Hcn4 + sinoatrial node (san, green), right atrium (pink, ra) and superior vena cava (gray, svc) from a WT control adult heart (left) in comparison to (B) a heart obtained from a littermate lacking both copies of an SAN-specific enhancer (ISE) for the transcription factor Isl1. Hearts were stained in whole mount with Hcn4 antibody and san, ra, and svc were manually segmented. Note the markedly reduced size of the Hcn4 + SAN in the Isl1 ISE/ ISE heart. (C) Electrocardiograms recorded from a WT control adult mouse (left) and (D) an Isl1 ISE/ ISE littermate demonstrate an example of sinus arrhythmias in animals lacking the Isl1 enhancer.
genome of mice to test whether target gene expression and SAN development are affected. In Galang et al., a 2.7-Kb segment at the Isl1 locus (hereafter, Isl1-Locus SAN Enhancer, or ISE) was deleted using CRISPR-Cas9. Homozygous Isl1 ISE/ ISE mice were born in Mendelian ratios and were viable and fertile with no major cardiac structural abnormality by echocardiography. However, detailed molecular and histological studies showed that Isl1 protein and mRNA expression were reduced specifically in Isl1 ISE/ ISE pacemaker cells, with a corresponding hypoplasia of the SAN and reduced pacemaker cell number (Figures 4A,B). Furthermore, heart rhythm monitoring revealed a slower heart rate in the mutant mice, as well as episodes of sinus arrhythmia and bradycardia, indicative of abnormal SAN function (Figures 4C,D). Taken together, these findings support a model in which ISE supports Isl1 transcription during SAN morphogenesis to ensure that the SAN achieves its normal size and cellularity, allowing for normal SAN function.
van Eif et al. deleted larger genomic regions within the SHOX2 and TBX3 TADs that contained differentially accessible ATACseq peaks. Examination of knockout embryos from a 250-kb deletion mouse line at the Shox2 locus showed loss of Shox2 expression in the SAN, with dysregulation of Shox2 target genes, hypoplasia of the SAN, and embryonic lethality presumed due to bradycardia. In a parallel experiment, a 280-kb putative regulatory region containing differentially accessible ATAC peaks was also deleted from the Tbx3 locus. Homozygous mutant mice had absent Tbx3 expression in the SAN and exhibited perinatal lethality, while heterozygous mice survived to adulthood and exhibited prolonged SAN recovery times with programmed stimulation, indicating an important role for this enhancer in regulating SAN electrophysiology.

Building a Model for Transcriptional Regulation in Pacemaker Cells Based on Loss of Function Studies and Epigenetic Profiling
By combining data from conditional loss of function studies for critical TFs, available ChIP-seq datasets for these factors, and the insights gleaned from the recent epigenetic profiling studies, a working model for transcriptional regulation in pacemaker cells can be constructed (Figure 5). Based on loss of function studies, Shox2 and Isl1 act upstream of Tbx3, which downregulates the atrial gene program. Shox2 also represses Nkx2.5 and its target genes. As noted above, within SAN pacemaker cell-specific ATAC-seq peaks, Isl1 binding motifs are enriched, suggesting it acts primarily as an activator of the SAN gene program. Outside of the SAN, Nkx-2.5 and Pitx2c repress the differentiation of additional pacemaker cells, further localizing the primary pacemaker of the heart specifically to the SAN. Overlapping ChIP-seq data for Gata4, Mef2c, and Tbx5 with SAN-specific ATAC-seq peaks confirms the binding of these TFs to putative SAN enhancers as well as more broadly to cardiac enhancers (Galang et al., 2020). In vivo studies have also demonstrated Mef2's upstream role in regulating both SAN and non-SAN genes, similar to Tbx5 and Gata4 (Vedantham et al., 2013). Remaining questions include the identities of factors that may Transcriptional control of the SAN gene program and right atrial cardiomyocyte gene program are contrasted. Both tissues require the cardiac factors Mef2, Gata4, and Tbx5 to active cardiac gene expression. However, reciprocal repression of SAN genes by Nkx2.5 and atrial genes by Tbx3, Tbx18, and Shox2, along with activation of SAN genes by Isl1 (with possible interacting partners) allows for cell type specific gene expression patterns. Question marks indicate a hypothesized interaction that has not been tested yet. Solid lines represent direct promotion/inhibition and dashed lines represent either direct or indirect promotion/inhibition. cooperate with Isl1 to regulate transcription in pacemaker cells. Since Isl1 is broadly expressed in other non-cardiac tissues, it is unlikely that Isl1 alone is sufficient to activate the pacemaker cell gene program. Candidate genes include Shox2 and Meis1, whose motifs are overrepresented near Isl1 sites in pacemaker cell-specific enhancers, as well as core cardiac TFs such as Tbx5, Gata4, Mef2, and others. In the coming years, as low-input ChIP-seq methodologies such as CUT&TAG are deployed more broadly, answers to these questions will be forthcoming.

RELATIONSHIP OF GENE REGULATION IN PACEMAKER CELLS TO HUMAN VARIATION AND DISEASE
Ultimately, a major goal of dissecting regulatory networks in pacemaker cells is to gain a better understanding of the pathophysiology of sinus node dysfunction and, more generally, to provide insight into the genetic underpinnings of variation in human heart rhythm. A number of previous studies have related genetic variation in heart rhythm parameters to regions of non-coding DNA that were subsequently found to contain enhancer that regulate expression of nearby genes (Christophersen et al., 2017;Roselli et al., 2018). Indeed, genome wide association studies (GWAS) have identified numerous single nucleotide polymorphisms (SNPs) in non-coding DNA that confer susceptibility to atrial fibrillation. Most notably, the chromosome 4q25 genomic locus, which exhibits the highest association with atrial fibrillation, contains enhancers that control expression of the homeodomain transcriptional factor Pitx2, a key regulator of left atrial identity and pulmonary venous (B) shows alignment of a resting heart rate GWAS at the TBX3 locus, with highly significant associations seen in the region with differentially accessible ATAC-seq peaks, underscoring the relevance of SAN enhancers to human sinoatrial node function. GWAS results were downloaded from http://www.nealelab.is/uk-biobank/. myocardium that also inhibits growth of the left-sided SAN node primordium (Roselli et al., 2018).
In keeping with these findings, GWAS that have been performed using resting heart rate or heart rate recovery after exercise have found numerous "hotspots" in non-coding regions of the genome that were presumed to have regulatory function (Kilpeläinen, 2016;Ramirez et al., 2018). One such region is closest to the gene MED13L, but also lies within a TAD that encompasses TBX3. SNPs in this region display strong associations with heart rate recovery after exercise. Not surprisingly, the GWAS hot spot overlaps all of the TBX3 enhancers that were identified and validated in van Eif et al. (2020), providing a conclusive demonstration of the relevance of epigenetic profiling to human heart rhythm (Figure 6). In the absence of the ATAC-seq data derived from pacemaker cells and the subsequent confirmatory studies, the association of a gene desert near MED13L with heart rate recovery would remain merely speculative. Galang et al. worked from the presumption that ATAC-seq peaks were more likely than other regions of the genome to harbor heart rate-associated SNPs, prompting them to perform a new association study on UK Biobank participants, limiting the analysis to SNPs occurring within 500 bp on either side of the human genomic regions syntenic to mouse pacemaker cell ATAC-seq peaks. The results of this analysis demonstrated many associations between heart rate and SNPs occurring near ATACseq peaks that are close to established genes that regulate heart rate. Furthermore, by focusing the association study on a subset of SNPs, and hence reducing the number of hypotheses tested, the uncorrected p value threshold needed to establish a significant association was reduced. Thus, by focusing only on differentially accessible peaks that are open in pacemaker cells, Galang et al. (2020) uncovered a significant association between heart rate and SNPs in the vicinity of the Isl1 SAN enhancer described above. Although the effect size was small, alongside the rest of the association study, this finding establishes that enhancers are important genomic loci for determination of sinus node function, and that epigenetic profiling techniques such as ATAC-seq have the potential to illuminate the genetic basis for variation in heart rhythm at the population level.

OUTSTANDING QUESTIONS AND FUTURE DIRECTIONS
The data summarized in this article have brought us closer than ever before to a detailed understanding of how pacemaker cells acquire their unique phenotypes during development, and to connecting the genetic determinants of pacemaker cell biology to vertebrate evolution and human physiology. Nevertheless, critical knowledge gaps remain.
First, while ATAC-seq data establish the genomic locations of the regulatory elements important for pacemaker cell-specific gene expression and indicate some of the transcription factors that are likely to bind to them, further work is needed to establish the precise timing and location of transcription factor binding events onto these enhancers. For example, which enhancers are specifically bound by Isl1 and what are its binding partners? What are the pioneer factors that open chromatin to allow for execution of the pacemaker cell gene expression program and what are the signaling pathways that activate these factors during sinus node development? Are these same pathways required for maintenance of pacemaker cell identify in the adult heart? How is epigenetic and transcriptional control related to age-related SAN dysfunction? Answers to these questions will fill in the missing pieces in our model for transcriptional control in pacemaker cells and will point the way toward translational strategies that could reverse SAN failure or even regenerate pacemaker cells.
A second area that will require further exploration is the connection between human genetic variation and SAN disease. As discussed, genome-wide and more limited association studies have demonstrated convincingly that a significant portion of human variation in SAN behavior is attributable to genetic differences at enhancer loci that regulate expression levels of critical SAN genes. A remaining challenge will be to take these analyses a step further by exploring whether variation in human non-coding regions is associated with the development of age-related or premature-onset sinus node dysfunction, atrial arrhythmias, and need for pacemaker implant. Defining the nature of these associations and the regulatory pathways involved could lead to new targets for interventions to prevent or treat SAN disease and other atrial arrhythmias.
Historically, research on these aspects of pacemaker cell and SAN biology has been challenging because of the difficulty of isolating pacemaker cells in large numbers, the lack of a suitable in vitro system to model pacemaker cells, and the limited access to human primary pacemaker cells. Fortunately, new tools and approaches are rapidly becoming available to surmount these technical hurdles. First, the development of human iPSC-based in vitro differentiation protocols, as highlighted in the study by van Eif et al., allows for the possibility of exploring the functions of cis-regulatory elements and transcriptional regulators in a more realistic human context without the limits imposed by low cell numbers. Thus, techniques such as ChIP-seq and Hi-C will eventually be used to create detailed, functionally annotated genomic maps of pacemaker cells. These datasets can then be integrated with data emerging from increasingly large cohorts of genotyped patients, including those with wholegenome sequencing alongside clinical data, to provide us with a detailed high-resolution model for how the epigenome and transcriptional biology connect with disease.
Alongside these technical innovations, the development of single cell genomic profiling technologies, including scRNA-seq and scATAC-seq, will circumvent the need to purify primary pacemaker cells for molecular analysis. Furthermore, data using these techniques can now be integrated across different platforms, allowing for more powerful in silico examination of cellular subtypes and cellular differentiation. Used in combination with new technologies that can isolate individual nuclei from differentiated tissues, these techniques will be invaluable as individual genes and pathways are explored in loss-and gain-of-function animal models and in small samples derived from human tissue. Taken together, the combination of recent technical and scientific advances leaves the field poised for rapid progress in the next several years.

CONCLUSION
Over the last several decades, an integrated model for the electrophysiological basis of pacemaker cell automaticity and SAN physiology has taken shape, but a detailed understanding of the underlying gene regulatory networks has remained elusive. Recent advances in developmental biology, chromatin biology, bioinformatics, and human genetics have now revealed that several key transcriptional regulators and genomic loci are critical for pacemaker cell development, SAN formation, and SAN function, and have clarified their relevance to human variation. Work in the coming years will translate these findings to an improved understanding of SAN disease, novel targets for therapeutics, and possibly SAN regeneration.

AUTHOR CONTRIBUTIONS
RM and VV wrote the manuscript and prepared the figures. CJ prepared the figures, performed the imaging and image processing, and edited the text for intellectual content. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
The figures were created with Biorender.com. 3D image reconstructions were created using the Nikon Imaging Center at UCSF.