In Vivo Dissection of Chamber-Selective Enhancers Reveals Estrogen-Related Receptor as a Regulator of Ventricular Cardiomyocyte Identity

Mammalian hearts comprise 2 major chamber types, atria and ventricles. The atria collect blood with high elasticity and propel it at relatively low pressure into the ventricles, which efficiently pump blood at high pressure.¹ To achieve these distinct functions, atrial and ventricular cardiomyocytes (aCMs and vCMs) have key morphological and functional differences that are encoded in distinct, chamber-specific gene expression programs.^2,3 One consequence is that atria and ventricles are primarily affected in distinct diseases. For instance, atria are directly affected by atrial fibrillation, whereas ventricles are more directly affected in dilated cardiomyopathy.^4–7 The transcriptional mechanisms responsible for chamber-specific cardiomyocyte gene expression remain largely unexplored, limiting our understanding of heart development and disease, and our ability to target genetic manipulations to aCMs or vCMs.

Mammalian transcriptional activation is achieved through the coordinated binding of transcription factors (TFs) to DNA cis-regulatory elements, including activating elements known as enhancers. TFs recruit other transcriptional regulators, including the transcriptional coactivator P300 into proximity with gene promoters through chromatin looping.^8–10 Cell specificity is achieved through regulation of each of these steps, including cell type–dependent expression of TFs and transcriptional regulators, enhancers accessibility, and enhancer looping.

Genome-wide association studies (GWAS) demonstrate significant association of genetic variation in noncoding regions to cardiovascular disease. Advancing from genetic association to mechanistic insight has been hindered by the challenges of functionally analyzing noncoding variation. In the context of heart disease, fully leveraging GWAS results would be facilitated by maps of noncoding regulatory elements active in aCMs and vCMs, including those selectively active in each of these cardiomyocyte subtypes.

Candidate transcriptional enhancers can be identified through their characteristic chromatin features, including acetylation of histone H3 at lysine 27 (H3K27ac),¹¹ P300 occupancy by P300,^12,13 which deposits H3K27ac, multiple TF co-occupancy,^14,15 and increased chromatin accessibility.^16,17 Although these high-throughput methods predict candidate enhancers, no feature or combination of features identifies enhancers with robust sensitivity and specificity. Moreover, these methods are not well suited to dissection of sequence features essential for enhancer activity. The massively parallel reporter assay (MPRA) has emerged as a high-throughput method to measure enhancer activity and to probe the sequences required for a region’s activity.^18,19 MPRA has been adapted to in vivo measurement of enhancer activity through the use of adeno-associated virus (AAV) to deliver MPRA libraries.^13,20

In this study, we mapped the chromatin features of aCMs and vCMs and nominated candidate chamber-selective enhancers on the basis of differential features. Candidate enhancers were tested in vivo using AAV-MPRA, leading to identification of 229 chamber-selective enhancers (CSEs). We characterized chromatin features of CSEs and used dense mutagenesis to identify their essential features. These studies suggested that estrogen-related receptor (ERR) promotes ventricular CSE activity. We validated this prediction by showing that ERR inactivation leads to the loss of vCM identity. Together these studies yielded a rich resource of chamber-selective chromatin features and CSEs and began to unravel the molecular basis for chamber-selective transcriptional programs.

Methods

Refer to the Supplemental Material for detailed Methods. Animal experiments were performed following protocols approved by the Boston Children’s Hospital Animal Care and Use Committee.

High-throughput data associated with this study, summarized in Extended Data 1, have been deposited in Gene Expression Omnibus (GSE215065). Data that support the findings of this study that are not within this repository or within the article are available from the corresponding author on reasonable request.

Statistical analyses were performed in Graphpad Prism, R, or python. Results are displayed as mean±SD. For box plots, the center line and hinges represent the 50th, 25th, and 75th percentile. The whiskers extend to the values within 1.5 times the interquartile range. Statistical tests used are specified in each figure legend. P<0.05 was considered statistically significant.

Results

Distinct Phenotypes and Transcriptional Programs of aCMs and vCMs

Differences in aCM and vCM form and function have been well described.^21–23 A comparison of the morphological features of aCMs and vCMs highlighted these differences. aCMs were shorter and thinner and had higher length-to-width ratios than vCMs (Figure S1A and S1B). This in situ imaging also showed that T tubules, tubular invaginations of the plasma membrane central to excitation-contraction coupling, are highly organized in vCMs and less prominent in aCMs.²³ In situ imaging with the mitochondrial dye tetramethylrhodamine, methyl ester, visualized the striking difference in mitochondria between aCMs and vCMs (Figure S1C).²²

We defined transcriptional differences that underlie these distinct phenotypes by performing RNA sequencing on purified neonatal aCMs and vCMs (Figure 1A; Extended Data 1 and 2). In total, 1126 and 872 genes had biased expression between aCMs and vCMs, respectively (P_adj<0.05 and |log₂ fold-change| > 0.58; Figure 1A). Gene ontology (GO) analysis of these genes indicated that aCM-biased genes are enriched in biological processes related to ion transmembrane transport, muscle contraction, and cell projection organization (Figure 1B), which may reflect the differences in geometry and T tubulation between aCMs and vCMs (Figure S1A and S1B). In contrast, the top 10 GO biological process terms for vCM-biased genes were related to metabolism and cellular respiration, and the top GO cellular component terms for vCM-biased genes were mitochondrial inner membrane and respiratory chain complex (Figure 1B), consistent with mitochondrial imaging (Figure S1C) and analyses of the human heart proteome.²⁴ Comparison with maturation-regulated ventricular genes (P28 versus P0) demonstrated that these regulatory programs contain some shared genes, especially enrichment of mitochondrial genes among both mature and ventricular cardiomyocytes, but are largely distinct (Figure S2).

• Download figure
• Download PowerPoint

These data identify transcriptional differences between aCMs and vCMs that contribute to functional differences in muscle contraction, ion transmembrane transport, and metabolism.

Cardiac TF and P300 Chromatin Occupancy in Atria and Ventricles

Cell type–specific gene expression is regulated by cell type–specific enhancers.⁸ We mapped atrial and ventricular enhancers by performing bioChIP-seq of transcriptional coactivator P300, a well-established marker of active enhancers,^12,13 and 7 key cardiac transcription factors (GATA4, MEF2A, MEF2C, NKX2-5, SRF, TBX5, TEAD1),^14,15 using mice in which FLAG and biotin acceptor peptide (BIO) epitopes were knocked into the C terminus of the endogenous genes (Figure 1C). Escherichia coli biotin ligase BirA, expressed from the Rosa26 locus, specifically biotinylated BIO, so that the labeled proteins could be pulled down on immobilized streptavidin with high sensitivity and specificity.^13,15,25 We acquired bioChIP-seq of these 8 transcriptional regulators in biological duplicate in neonatal atria and ventricles (Figure 1D; Extended Data 1 and 3). High data quality was indicated by high correlation of biological replicates (Figure S3) and by motif enrichment analysis (Figure 1E and Figure S4A). Each TF bioChIP-seq enriched for biological process GO terms that differed between TFs and between chambers (Figure S4C).

The 7 TFs yielded 152 832 reproducible peaks, of which 27 930 and 60 465 were called only in atria or ventricles, respectively (Figure 1D and Figure S5A; Extended Data 3). P300 bioChIP-seq yielded 13 386 reproducible peaks in atria (aP300) and 29 804 in ventricles (vP300), with 1145 and 8368 peaks that were specific to atria or ventricles, respectively (aSpP300 and vSpP300; Figure 1D and Figure S5B; Extended Data 3). We previously showed that regions co-occupied by multiple TFs (≥5 TFs) are likely to be active enhancers, and that these regions only partially overlap with P300 regions.^14,15 By integrating the TF bioChIP-seq data, we identified 2423 and 4196 regions occupied by multiple TFs (MTFs: ≥5 TFs) in atria (aMTF) or ventricles (vMTF), respectively (Extended Data 4). Thirty-five percent and 18% of aMTF and vMTF regions, respectively, lacked P300 co-occupancy (Figure 1F), consistent with our previous studies.^14,15 There were 454 and 1666 MTF regions specific to atria or ventricles (aSpMTF or vSpMTF, defined as bound by ≥5 TFs in 1 chamber and ≤3 TFs in the other chamber; Figure S5A and Extended Data 4).

Taken together, bioChIP-seq reproducibly mapped chromatin occupancy by P300 and cardiac TFs in atria and ventricles. These data identified chamber-specific P300 or MTF regions that are candidate CSEs.

CSEs Identified in Vivo by AAV-MPRA

To measure the activity of candidate CSEs in aCMs and vCMs, we deployed an AAV-based MPRA¹⁵ on the basis of the self-transcribing active regulatory region sequencing design,¹⁹ in which an enhancer is placed in the 3′-UTR of an mCherry reporter gene, downstream of an hsp68 minimal promoter. We confirmed that this vector could detect CSE activity by individually testing 25 regions with chamber-selective TFs or P300 binding using a modified vector that also contained an RNA Polymerase III (U6) promoter driving a reporter RNA, Broccoli, to facilitate normalization for AAV transduction efficiency (Figure 2A). After systemic delivery, we measured enhancer activity by mCherry imaging. Of 25 tested candidates, 13 drove detectable mCherry expression and 9 exhibited chamber-specific activity (Figure 2B and 2C and Figure S6). Atrium-selective enhancers neighbored Myl7, Fgf12, Gpx3, and Bmp10, and ventricle-selective enhancers neighbored Myl3, Myl2 (2 separate regions), Irx1, and Irx4 (Figure 2B and Figure S6A and S6B). Chamber selectivity of 8 of these 9 enhancers was validated by reverse transcriptase-quantitative polymerase chain reaction for mCherry, normalized to Broccoli (Figure 2C and Figure S6C and S6D).

• Download figure
• Download PowerPoint

We next selected candidate regions to screen for CSE activity. Our previous data showed that regions bound by different numbers of TFs had different sensitivity and specificity to predict enhancers. We nominated 2943 candidate regions belonging to the following 4 categories (Figure 2D, Figure S5, and Extended Data 5): (1) Regions bound by 1 to 2 TFs (TF1-2); (2) regions bound by 3 to 4 TFs (TF3-4); (3) regions bound by ≥5 TFs (MTF); and (4) regions bound by P300. For MTF and P300 regions, we defined a subset with chamber-selective occupancy (aSpMTF or aSpP300; vSpMTF or vSpP300). We also included 954 negative control regions (P300 regions in murine embryonic stem cells but not heart¹³). Using pooled oligonucleotide synthesis, each 400-bp region was synthesized as 230 nucleotides (nt) self-priming oligonucleotide pairs (Figure 2D). The region pool was cloned into the AAV-MPRA and packaged into an AAV library. The pooled AAV library was delivered to P0 pups. Atria and ventricles were harvested at P7 with 5 atrial and 5 ventricular replicate samples. The frequency of each region was determined in reverse-transcribed RNA and AAV genome DNA. Candidate regions with low coverage in AAV genome DNA were removed (Figure S7A). For the surviving 2160 candidate and 479 negative control regions, enhancer activity was quantified by the ratio of reads from RNA compared with DNA. The replicates from the same chamber clustered together with high correlation, indicating high reproducibility of the assay (Figure S7B).

In total, 1092 enhancers had activity in the AAV-MPRA, where active regions were defined as those with significantly increased representation in RNA compared with DNA (Figure 2E and 2F; Extended Data 5). Some 966 regions were active in atria and 963 in ventricles. Regions bound by more TFs were more likely to be active enhancers in each chamber. aMTF/vMTF or aP300/vP300 regions were most enriched for MPRA activity, whereas negative control regions were depleted for active regions. The level of enhancer activity in atria and ventricles was greatest for aMTF/vMTF and aP300/vP300 categories. TF1-2 did not significantly differ from negative control regions, and TF3-4 was slightly greater than negative controls in atria but not ventricles (Figure 2E and 2F).

Analysis of active enhancers for chamber-biased activity (|log₂ A_activity/V_activity|>0.58 and false discovery rate<0.05). identified 110 atrium-specific enhancers (aSpEnh) and 119 ventricle-specific enhancers (vSpEnh; Figure 2G). The aSpEnh were enriched for aSpMTF and aSpP300 regions, and vSpEnh were enriched for vSpMTF and vSpP300 regions, whereas regions bound by MTF or P300 in both chambers (non-sel) were not enriched among either aSpEnh or vSpEnh (Figure 2G). Both aSpMTF/vSpMTF and aSpP300/vSpP300 regions had significantly greater chamber selectivity than the regions without selective occupancy. These data indicate that candidate CSEs can be identified from regions with chamber-specific MTF or P300 occupancy, although only a subset of such regions had chamber-selective activity.

GWAS have implicated noncoding regions in the modulation of human cardiac phenotypes. We integrated 521 unique lead single-nucleotide polymorphisms with 10 866 corresponding LD single-nucleotide polymorphisms (r²>0.9, within 50 kb of lead single-nucleotide polymorphisms) from cardiovascular GWAS studies to define loci related with cardiac phenotypes or diseases (Table S1). Sixteen of these loci contained 18 MPRA-defined active enhancers, of which 4 were chamber selective (Table S2). Twelve, including loci near TBX5, MICU2, ITGA9, and ATP2A2, were associated with ECG features, 1 near CAV1 with atrial fibrillation, 1 near C6orf15 with dilated cardiomyopathy, and 1 near PLPP3 with myocardial infarction (Figure 2E and 2F). These findings suggest a potential link between these enhancers and human cardiovascular traits and diseases.

Enhancers recruit TFs to regulate spatiotemporal gene expression. We identified motifs enriched in active enhancers to implicate regulatory TFs. Distinct sets of motifs were enriched among aSpEnh and vSpEnh (Figure 2H). Compared with genomic background, the core cardiac TF motifs GATA4, MEF2, TEAD1, and NKX2-5 were enriched among both groups of CSEs and the T-Box motif was enriched only among aSpEnh. To better pinpoint motifs responsible for chamber-selective activity, we used CSEs from the opposite chamber as the background for motif enrichment analysis. TBX5, MEIS1, and PITX2 motifs were among those enriched in aSpEnh compared with vSpEnh. These TFs regulate atrial ion channel genes and are implicated in atrial arrhythmias and atrial restrictive cardiomyopathy.^26–29 In contrast, vSpEnh were most highly enriched compared with aSpEnh for motifs of MEF2, SF1, and ERR and other nuclear receptors (EAR2, NR5A2). ERR stimulates expression of genes involved in mitochondrial biogenesis and oxidative phosphorylation,^30,31 providing a link between ventricular CSEs and metabolic genes enriched among vCM-biased genes (Figure 1B).

Taken together, AAV-MPRA testing of 2639 candidate cardiac enhancers identified 1092 active regions, including. 229 that were atrial or ventricular CSEs. CSEs were enriched for chamber-specific MTF and P300 occupancy. Atrial and ventricular CSEs were enriched for TBX5 and ERR motifs, respectively.

Distinct Chromatin Accessibility of CSEs Between aCMs and vCMs

Cell type–specific chromatin accessibility is involved in determining and maintaining cellular identity.^17,32 We measured chromatin accessibility of purified neonatal aCMs and vCMs in biological duplicate by ATAC-seq. There was high correlation between the replicate samples (Figure S8A; Extended Data 1 and 3); 59 939 accessible regions were reproducible across the 2 chambers; 4060 were atrium specific and 12 970 were ventricle specific (Figure 3A).

• Download figure
• Download PowerPoint

We then analyzed the chromatin accessibility of CSEs identified by AAV-MPRA. aSpEnh and vSpEnh were clearly separated into 2 groups on the basis of their ATAC signal in aCMs and vCMs (Figure 3B). The chromatin accessibility of aSpEnh regions was significantly higher in aCMs than vCMs, and, reciprocally, accessibility of vSpEnh regions was significantly higher in vCMs than aCMs. In contrast, accessibility of regions without detectable enhancer activity or regions with enhancer activity in both aCMs and vCMs showed no significant difference between aCMs and vCMs (Figure 3C). For example, an aSpEnh adjacent to atrial gene Myl7 showed atrium-selective occupancy by P300 and MTFs and had strong, aCM-selective accessibility. In contrast, a vSpEnh adjacent to ventricular gene Myl3 had ventricular-selective occupancy by P300 and MTFs and strong, vCM-selective accessibility (Figure 3D).

Taken together, our data indicate that chamber-selective chromatin accessibility plays an important role in regulating CSE activity and selectivity.

Chamber-Selective Chromatin Loops Link a Subset of CSEs to Target Genes

Enhancers can regulate genes several hundred megabytes away, with chromatin loops positioning distal regulatory elements near their target promoters.^10,33,34 To measure the contribution of 3-dimensional genome structure to chamber-specific transcriptional programs and CSE activity, and to identify target genes linked to CSEs, we used H3K27ac HiChIP³⁵ to identify enhancer chromatin contacts in purified neonatal aCMs and vCMs (Extended Data 1 and 6). Biological triplicate replicates grouped together, with clear separation between aCMs and vCMs (Figure S8). We identified 80 792 high-confidence, reproducible chromatin loops in aCMs and 53,840 loops in vCMs (Figure 4A), with 3457 (4.3%) and 2150 (4.0%) involving a region with atrial or ventricular MPRA activity, respectively (Figure 4A).

• Download figure
• Download PowerPoint

We visualized H3K27ac HiChIP loops involving CSEs near atrial gene Myl7 and ventricular gene Myl2 (Figure 4B). We observed higher interaction frequency in aCMs near Myl7 and in vCMs near Myl2. Genome browser views at these genes showed that the Myl7 aSpEnh had robust atrium-specific H3K27ac signal and strong chromatin interaction with the Myl7 promoter (Figure 4C, top). Similarly, the Myl2 vSpEnh exhibited high ventricle-selective H3K27ac signal and strong chromatin interaction with the Myl2 promoter (Figure 4C, bottom).

In total, there were 4643 (5.7%) and 5420 (10.1%) loops specific to aCMs or vCMs, respectively. MPRA-identified chamber-selective enhancers were enriched for overlap with chamber-selective loops (atria odds ratio=4.04, Fisher P=1.2E-37; ventricles odds ratio=2.76, Fisher P=3.9E-12; Figure 4D). However, most MPRA-identified chamber-selective enhancers (atria, 63/110; ventricles, 78/119) did not overlap chamber-selective loops. These data indicate that chamber-selective looping is enriched at CSEs but is not a predominant mechanism for CSE regulation.

Using the H3K27ac HiChIP atlas to associate enhancers to promoters in aCMs and vCMs, we assessed the contribution of CSEs to chamber-selective transcriptional programs. Anchor genes (genes with promoters within loop anchors) of chromatin loops in aCMs and chromatin loops in vCMs each had slightly biased expression toward aCMs and vCMs, respectively (Figure 4E). Anchor genes of specific chromatin loops in aCMs had significantly more biased expression toward aCMs, compared with anchor genes of chromatin loops in aCMs. vSpLoop anchor genes also had vCM-biased expression, although this did not reach statistical significance. This analysis indicates that chamber-selective promoter looping contributes to chamber-selective gene expression.

We calculated the sensitivity and specificity of MTF, P300, open chromatin, and chromatin looping for identifying MPRA active regions (Figure 4F). MTF occupancy had the greatest specificity (≈75%), and P300, open chromatin, and looping had moderate specificity (42%–63%). Each of the 4 features also had moderate sensitivity (37%–68%). The union of the 4 parameters had high sensitivity (89%–93%), indicating that most active enhancers have at least 1 of these features. A parallel analysis focusing on the ability of chamber-specific features to discern CSEs from active enhancers showed that chamber-specific MTF, P300, and open chromatin had comparably high specificity (>93%), whereas chamber-specific looping was slightly less specific (≈80%; Figure 4G). However, each feature and the union of features had low sensitivity (<26%) for detecting CSEs among active enhancers, suggesting the importance of other features (eg, TF motifs) to determine chamber-specific activity.

Systematic Mutagenesis to Identify Features Required for CSEs

To investigate the sequence features required for CSE activity and selectivity, we performed dense mutagenesis of 13 atrium-selective, 8 ventricle-selective, and 8 shared enhancers. Each of these 29 400-bp candidate regions was synthesized as 80 180-nt oligonucleotides tiled in 5-nt steps (wild-type, WT; Figure 5A). Each WT 180-nt oligonucleotide was mirrored by a counterpart in which the central 5 nt was deleted (Mut; Figure 5A). Each sequence was identified by a unique 10-nt barcode. From the initial MPRA library, we included negative control regions, inactive in both chambers, and 135 positive control regions, active in both chambers, yielding a MPRA library of 5274 oligos (Extended Data 7). We constructed a pooled AAV library. After excluding oligos with low coverage (<5 fragments per million; Figure S9A), the library contained 2394 WT sequences and 1761 WT-Mut pairs. The library was delivered to P0 pups and enhancer activity of each sequence was analyzed at P7 by amplicon sequencing. Enhancer activities were well correlated across 5 atrial and 4 ventricular replicates (Figure S9B).

• Download figure
• Download PowerPoint

We analyzed the activity of WT 190-bp “tiles” for their activity in aCMs and vCMs to identify the subset of original 400-bp enhancers whose activity was reproduced in 180-bp subregions (Figure 5B, Figure S10A and S10B, and Extended Data 1 and 7). Active tiles clustered together (Figure 5B, red areas), highlighting the “active core” of the original 400-bp enhancers that were sufficient for activity. We found that 11 of 13 aSpEnh, 7 of 8 shared, and 7 of 8 vSpEnh retained activity in at least 3 consecutive 180-bp WT tiles.

Next we compared the atrial and ventricular enhancer activity of each WT-Mut pair (Figure 5C and Figure S10C and S10D). This identified 5-bp “bins” that were required for region activity or selectivity. These tended to localize within the active core of each active enhancer. A detailed view of the dense mutagenesis of the ventricle-selective Myl3 enhancer is shown as an example (Figure 5D). The Activity track identified the left side of the enhancer as its active core, which exhibits strong ventricular and weak but detectable atrial enhancer activity. Mutant tiles within this core identified several bins that are essential for both levels of activity and chamber selectivity (arrowheads, Figure 5D). Bins 18 to 19 contain a MEF2 motif whose deletion strongly reduced enhancer activity in both atria and ventricles (tracks Mut/WT A and V) and strongly reduced ventricular selectivity (tracks V/A-Mut versus WT). Bin 12 contains an ERR motif whose deletion likewise strongly reduced enhancer activity in both atria and ventricles and strongly reduced ventricular selectivity. At the edges of the active core, bin 26 containing an EBF1 motif also impacted enhancer activity in both aCMs and vCMs. Together, the systematic mutagenesis identified MEF2 and ERR binding sites that are essential for the ventricle-selective activity of this enhancer.

We performed a similar analysis of the systematic mutagenesis of other atrium- and ventricle-selective enhancers (Figure S10E–S10G). For example, an atrial-selective enhancer (chr6:30 545 609–30 546 008) linked by H3K27ac HiChIP to Mir335, an atrial fibrillation–associated microRNA,³⁶ required the T-box motif in bins 43 to 44 and 51 to 52 (Figure S10F) for atrial selectivity (V/A-Mut track), and deletion of bin 50 significantly increased enhancer activity in ventricle (Mut/WT A and V tracks). Thus, systematic mutagenesis of this enhancer identified T-box binding sites essential to its atrial-selective activity, through enhancer activation in aCMs and repression in vCMs.

To identify TF motifs important across the mutagenesis experiment, we scanned mutant and WT pairs with significantly different activities for enriched motifs (see Methods). Pairs without significant change in activity were used as controls. Because a deletion can remove or insert a motif and increase or decrease enhancer activity, we considered 4 classes of mutations. Among mutations that caused loss of motif and loss of activity, the ERR motif was among the most highly enriched in ventricles (Figure 5E), consistent with the enrichment of the ERR motif among 400-bp enhancers with ventricle-selective activity (Figure 2H). Among mutations that caused loss of motif and increased activity, the TBX5 motif was the most enriched in ventricles (Figure 5E), suggesting that this motif contributes to chamber-selective enhancer repression in ventricles. In atria, however, mutations that created the TBX5 motif increased enhancer activity, suggesting that overall TBX5 acts as an activator in aCMs but a repressor in vCMs.

Taken together, dissection of chamber selective by systematic mutagenesis identified TF motifs that are essential for enhancer activity and selectivity in aCMs and vCMs. Most notably, the ERR motif was essential for ventricular enhancer activity.

ERRα/γ Regulate Ventricular Cardiomyocyte Identity

The 2 MPRA experiments indicated that the ERR motif is essential for ventricular enhancer activity (Figures 2H and 5E), including Myl3-Enh (Figure 5D). The ERR motif is bound by ERR family TFs, with ERRα and ERRγ being the predominant cardiac isoforms.^30,37 ERRα/γ were more highly expressed in vCMs than aCMs (Figure S11A and S11B), To directly test ERRα/γ requirement for ventricular enhancer activity, we treated ERRα^flox/flox; ERRγ^flox/flox mice with the AAV Myl3 enhancer-mCherry reporter construct (Figure 3G) and either AAV-Tnnt2-Cre or AAV-Tnnt2-GFP (Figure S11C and S11D). Double cardiac knockout (dCKO) mice had reduced ventricular Myl3 enhancer activity as determined by mCherry imaging and reverse transcriptase-quantitative polymerase chain reaction (Figure 6A and 6B). Moreover, endogenous Myl3 gene expression was also reduced in dCKO (Figure 6C). We also examined the effect of dCKO on the activity and selectivity of the Myl7 atrium-specific enhancer. This enhancer became ectopically active in dCKO ventricles (Figure 6D). reverse transcriptase-quantitative polymerase chain reaction confirmed 17.4-fold increased Myl7 enhancer activity in dCKO ventricles (Figure 6E). Accordingly, we also detected elevated ventricular expression of endogenous Myl7 (Figure 6F). These data suggested that ERR inactivation led to impaired ventricular identity, with reduced ventricular enhancer activity and increased atrial enhancer activity.

• Download figure
• Download PowerPoint

We assessed the effect of ERR inactivation on chamber-selective gene expression across the transcriptome by analyzing RNA-sequencing data from P35 ERRα/γ dCKO and control vCMs.³⁰ Using principal component analysis, we compared the effect of ERRα/γ dCKO with the difference in gene expression between P0 vCMs and aCMs (Figure 6G). The 2 datasets were separated along PC1, reflecting differences in both stage and technique between these datasets. PC2 separated vCMs from aCMs, and also control from dCKO, and indicated that dCKO cardiomyocytes are more similar to aCMs whereas control cardiomyocytes are more similar to vCMs.

We further analyze the effect of ERR inactivation on chamber-selective gene expression by performing Gene Set Enrichment Analysis³⁸ on the P35 ERRα/γ dCKO and control vCM RNA-sequencing dataset. vCM-biased genes were highly enriched among genes downregulated by ERRα/γ inactivation, and aCM-biased genes were highly enriched among upregulated genes (Figure 6H). Of the 422 vCM-biased genes that were differentially expressed in dCKO vCMs, 325 (77%) were downregulated, including Myl2, Myl3, Irx1, and Myh7b (Figure 6I, Fisher P=4.3E-36). On the other hand, of the 420 aCM-biased genes that were differentially expressed in dCKO vCMs, 298 (71.0%) were upregulated, including Sln, Fgf12, Myl7, Bmp10, Pitx2, and Nppa (Figure 6F, Fisher P=1.8E-18).

ERRα/γ regulate mitochondrial activity and biogenesis,³⁷ and vCM-biased genes are enriched for mitochondrial functional terms (Figure 1B). Of 1140 mitochondrial genes,³⁹ 577 (49%) were differentially expressed in ERRα/γ dCKO vCMs. After removing mitochondrial genes, the remaining vCMs-biased genes were still highly enriched among dCKO-downregulated genes (Figure S11E and S11F, Fisher P=8.3E-6). Enriched GO terms included muscle contraction, ion transport, and second messenger signaling (Figure S11G). These results indicate that ERRα/γ regulate elements of ventricular identity beyond mitochondrial biogenesis.

ERRα/γ regulate cardiomyocyte maturation.^30,40 Forty-six percent of genes dysregulated in dCKO vCMs overlapped with vCM genes with altered expression during postnatal maturation (Figure S12A),⁴¹ although the maturational expression change correlated poorly with the expression change in dCKO (R²=0.2437; Figure S12B). The mutually dysregulated genes were enriched for mitochondrial and metabolic GO terms, and muscle system process as well (Figure S12C). Moreover, after removing maturation differentially expressed genes, aCM-biased and vCM-biased genes were still enriched among dCKO-upregulated and -downregulated genes, respectively (Figure S12E–S12G). These results suggest that ERRα/γ regulates ventricular identity beyond its regulation of cardiomyocyte maturation.

Together, these data strongly confirm that ERR promotes ventricular identity.

Discussion

We investigated the transcriptional mechanisms that maintain chamber-selective gene expression in cardiomyocytes. Our atlas of chromatin occupancy by cardiac TFs and P300, accessible chromatin, and chromatin looping will be a rich resource for future studies of chamber-specific gene regulation. Using AAV-MPRA, we identified active and chamber-selective enhancers. Systematic mutagenesis identified motifs of key TFs required for activity of these CSEs; among these were TBX5 in aCMs and ERR in vCMs. Our studies of ERRα/γ dCKO vCMs validated the key role of ERR in maintaining ventricular identity, by promoting expression of ventricular genes and repressing atrial genes.

We studied the chromatin features of enhancers active in atria or ventricles. Enhancer features were enriched at MPRA-active regions and at CSEs. For identifying CSEs from active enhancers without chamber selectivity, chamber-selective MTF, P300, and ATAC were all highly specific, and chamber-selective enhancer loops were slightly less specific. However, none of the chromatin features or even their union was able to sensitively identify CSEs from non–chamber-selective active enhancers. This indicates that other features, such as individual TF motifs or combinations of motifs, are important determinants of CSEs.

We focused on factors that maintain chamber-selective gene expression rather than on atrial and ventricular specification and differentiation, the focus of prior studies. NKX2-5, MEF2C, and IRX4 have been implicated in ventricular cardiomyocyte differentiation.^42–44 We found that MEF2 was also required to maintain ventricular enhancer activity and selectivity. We did not detect NKX2-5 or IRX4 motifs, suggesting that the chamber-selective function of these TFs is restricted to developing cardiomyocytes. The most highly enriched ventricular enhancer motif belonged to ERR, and systematic tiling mutagenesis and transcriptome-wide analysis of ERRα/γ double-knockout vCMs confirmed that it is essential for promoting ventricular identity. This requirement extended beyond its previously reported functions to promote mitochondrial biogenesis or cardiomyocyte maturation.^30,37,40

TBX5 and NR2F2 are critical for atrial differentiation.^45,46 Ablation of NR2F2 later in embryonic development did not affect atrial identity,⁴⁶ suggesting that it is not required for maintenance of atrial identity. TBX5, on the other hand, remains essential for expression of many genes in mature atria, and adult-stage TBX5 deficiency causes atrial fibrillation.⁴⁷ Our analysis of atrial-specific enhancers supports the hypothesis TBX5 is essential for maintenance of atrial gene expression.

GWAS studies link noncoding sequence variants to cardiac phenotypes, likely through their effect on CRE activity. Eighteen of the murine aCM or vCM enhancers that we identified were contained within noncoding regions linked to cardiac diseases. For instance, active enhancers of Tbx5 and the key Ca²⁺ transporter Atp2a2 map to human loci associated with the electrocardiographic PR interval and QT intervals, respectively. The combination of in vivo MPRA assays, chromatin feature, and 3D genome maps will contribute to the ongoing functional dissection of clinically relevant regulatory loci and associated sequence variants.

Regulatory elements with well characterized properties are important reagents for research and therapeutic applications. The enhancers that we have identified, already tested in AAV, have a range of enhancer strength and chamber selectivity, providing a rich resource for the selection of elements with the desired regulatory properties. For example, the 110 aCM-selective and 119 vCM-selective enhancers that we identified could be used for chamber-selective expression in research studies or for therapeutic gene delivery.

Our study had some limitations. Our study addresses maintenance chamber–selective expression, rather than the specification or differentiation of aCMs or vCMs. We used neonatal cardiomyocytes for the chromatin feature mapping and MPRA assays, reasoning that most aCM and vCM differences had been established by this stage. This stage yields the highest quality chromatin data and is most amenable to functional studies. Our 3-dimensional genome data were based on a single method and a single antibody. Additional experiments are needed to validate our findings in adult cardiomyocytes, in homeostasis and in disease states.

Article Information

Acknowledgments

The authors thank Dr Johan Auwerx for sharing ERRα/γ flox mice. Drs Cao, X. Zhang, and Akerberg contributed equally to this work. Drs Cao, X/ Zhang, Akerberg, and Pu conceived and designed the study and analyzed data. Dr Akerberg performed bioChIP-sequencing, ATAC-sequencing, and RNA-sequencing experiments. Dr Cao performed the MPRA, H3K27ac HiChIP, and individual enhancer assays. Dr X. Zhang developed custom MPRA design and analysis pipelines. Drs Zhang and Cao conducted bioinformatic analyses. Drs Yuan, Sakamoto, Xiao, Van Dusen, Shou, Sweat, Y. Wang, Chen, Y. Zhang, and P. Wang contributed reagents, data, and analyses. Drs Cao, X. Zhang, and Pu wrote the manuscript with contributions from the other authors. Drs Pu and Kelly oversaw the project. All authors read and approved the manuscript.

Supplemental Material

Extended Methods

Figures S1–S12

Tables S1–S5

Extended Data 1–7

References 48–75

Footnotes