Tbx5 maintains atrial identity in post-natal cardiomyocytes by regulating an atrial-specific enhancer network

Abstract

Understanding how the atrial and ventricular heart chambers maintain distinct identities is a prerequisite for treating chamber-specific diseases. Here, we selectively knocked out (KO) the transcription factor Tbx5 in the atrial working myocardium to evaluate its requirement for atrial identity. Atrial Tbx5 inactivation downregulated atrial cardiomyocyte (aCM) selective gene expression. Using concurrent single nucleus transcriptome and open chromatin profiling, genomic accessibility differences were identified between control and Tbx5 KO aCMs, revealing that 69% of the control-enriched ATAC regions were bound by TBX5. Genes associated with these regions were downregulated in KO aCMs, suggesting they function as TBX5-dependent enhancers. Comparing enhancer chromatin looping using H3K27ac HiChIP identified 510 chromatin loops sensitive to TBX5 dosage, and 74.8% of control-enriched loops contained anchors in control-enriched ATAC regions. Together, these data demonstrate TBX5 maintains the atrial gene expression program by binding to and preserving the tissue-specific chromatin architecture of atrial enhancers.

Extended Data Figure 1.. Characterization of the cardiac expression domain of AAV9:Nppa-EGFP.

a-b, Hcn4^CreERT2/+; Rosa26^LSL-Tomato pups were injected with 2 X 10¹¹ viral genomes per gram bodyweight (VG/g) AAV9:Nppa-EGFP at P8 and injected with tamoxifen at P21 and P22. Hearts were harvested in PBS and brightfield, GFP, and Tomato fluorescent images were acquired. A similar result was observed with administration at P2. c, Hearts were sectioned and stained with WGA. Tomato signal marked the sinoatrial node (SAN) and atrioventricular node (AVN) but not the working myocardium. GFP was restricted to the working myocardium. A similar staining pattern was observed in 4 other littermates.

Extended Data Figure 2.. Phenotypic characterization of Tbx5^AKO mice.

a, Trichrome staining of Tbx5^Flox/Flox and Tbx5^AKO hearts. b, RTqPCR quantification of Postn (periostin) and Col1a1 (collagen type 1 alpha 1 chain), two indicators of fibrosis (n = 3 control and 4 KO mice). Unpaired 2-sided t-test. Error bars represent mean values +/-SEM. c-d, Immunostaining for sarcomeric α-actinin (SAA) or FSD2, markers of the Z line and the junctional Sarcoplasmic reticulum, respectively, in the left atrium of the indicated genotypes. Localization of both proteins is disrupted in Tbx5^AKO atria. e, Pattern of SAA signal intensity. SAA intensity along the long axis of cardiomyocytes demonstrated a periodic signal in control, consistent with regular position of sarcomere Z-lines, and loss of periodicity in Tbx5^AKO. f, Preserved ventricular function of Tbx5^AKO mice. Tbx5^Flox/Flox mice were treated with AAV9:Nppa-EGFP (control) or AAV:Nppa-Cre (Tbx5^AKO) at P2. Echocardiography was performed at P20. LVID;d, left ventricular internal diameter at end diastole. EF, ejection fraction. n = 5 control and 4 KO mice. Error bars represent mean values +/-SEM. g, Surface EKG recordings. Orange and black bars highlight successive RR intervals. h, Poincaré plots. The RR interval of greater than 1500 beats on EKG recordings is plotted versus the RR interval of the subsequent beat (RR[+1]). to visualize the dispersion of interbeat intervals in Tbx5^AKO, consistent with atrial fibrillation. i, Standard deviation of the RR interval, a measure of heart rate irregularity, was calculated for at least 1500 beats for each group at P21. n = 4 control and 9 KO mice. Unpaired two-sided t-test: ***, P=0.005. Error bars represent mean values +/-SEM. j, Time course of heart rate irregularity. Serial EKGs were acquired from control or Tbx5^AKO mice at the indicated time points. SDRR was measured and compared between groups using a two-way ANOVA. Sidek’s multiple comparison test was used to compare between genotypes at each time point. **, P<0.01. ***, P<0.001. For P14 and P17 timepoints, n = 3 control and 3 KO mice. For P8 and P21 timepoints, n = 3 control and 8 KO mice. Error bars represent mean values +/-SEM. k, Simultaneous intracardiac and surface ECG recordings demonstrate normal synchronous atrial-ventricular rhythm in TBX5Flox/WT mice injected with AAV9-Nppa-Cre. l, In contrast, animals with complete atrial ablation of TBX5 (AAV-Nppa-Cre + TBX5^flox/flox) demonstrate nearly continuous low-amplitude atrial activity and and irregular ventricular response consistent with atrial fibrillation.

Extended Data Figure 3.. Tbx5 overexpression atrializes ventricular myocytes.

a, Strategy to generate TBX5-OE ventricles. A full length Tbx5 cDNA downstream of the ubiquitous CAG promoter is activated by the cardiomyocyte specific Myh6-Cre transgene. b-c, Volcano plot comparing the change in gene expression of mouse left ventricle overexpressing Tbx5 compared to control LV. aCM genes are marked in red in (a) and vCM genes are denoted in blue in (b). b-c, Wald’s test followed by Benjamini Hochberg correction. d-e, Fisher’s exact test (two tailed) was performed to determine if changes in chamber selective gene expression downstream of Tbx5 overexpression were significant. f-h, Staining for MYL4, MYL7, and MYL2 in atria and ventricles of control and TBX5^OE hearts. Staining for each marker was performed using 3 control hearts and 2 TBX5 OE hearts, and representative images are shown.

Extended Data Figure 4.. Single cell dataset metrics.

a, Parameters from single nucleus datasets. Mean values per nucleus: nCount_ATAC, number of ATAC fragments; nFeature_ATAC, number of ATAC peaks with at least one read; nCount_RNA, number of RNA fragments; and nFeature_RNA, number of genes. b, Transcription start site (TSS) aggregation plots for the scATAC multiome datasets showing the expected enrichment of ATAC fragments. c, WNN UMAP of the dataset split by original sample. Each of the KO and control replicates have a high degree of overlap, demonstrating high reproducibility of cell state changes in Tbx5^AKO atria.

Extended Data Figure 5.. Differentially expressed genes between myocyte clusters.

a, Volcano plot of differentially expressed genes (DEGs) between early and late pseudotime clusters in the control trajectory (Myocyte_4 vs. Myocyte_1, left) and the KO trajectory (Myocyte_5 vs. Myocyte_6; right). Wilcoxon rank sum test , Bonferroni correction. b. MA plot of RNA-seq experiment comparing P0 and P28 aCMs. c-d. Volcano plots shown in (a) overlayed with genes enriched in P28 and P0 aCMs. Wilcoxon rank sum test , Bonferroni correction. e. The proportion of DEGs from comparisons in (a) that overlap genes selectively expressed in aCMs at P0 and P28. We did not observe enrichment of P0 or P28 selective aCM genes in early or late pseudotime clusters, respectively. This suggested that pseudotime trajectories did not correspond to chronological time. f-g. GO biological process terms enriched for DEGs for the comparisons shown in (a). The top 10 terms enriched by genes upregulated in the indicated cluster are shown. Functional terms related to cardiomyocyte function were enriched in the late pseudotime clusters, Myo_1 (control) and Myo_6 (AKO). Fisher’s Exact Test, Bonferroni correction. h, WNN UMAP plot colored by a “functional cardiac gene” index, which was calculated based on the aggregate expression of the six indicated genes, which are required for the efficient pumping function of aCMs.

Extended Data Figure 6.. A comparison of TBX5 RNA-seq datasets.

We compared snRNAseq data from Tbx5^AKO aCMs to three other bulk RNA-seq datasets involving Tbx5 gain- or loss- of function. Control_aCM indicates Myocyte_1 and KO_aCM indicates Myocyte_6. a, c, e: Scatter plot of fold-change in Tbx5^AKO aCMs compared to the other three datasets. Left, points are colored by significance in each dataset (P_adj<0.05). Right, genes with significant differential expression in both datasets are colored by chamber selectivity. Pie charts summarize the proportion of each class of genes within each quadrant. b, d, f: Enrichment of aCM or vCM genes in the indicated quadrants of the scatter plots. Fisher’s exact test, two-tailed. a-b. Comparison to LA tissue with ubiquitous, adult-induced inactivation of Tbx5 (Tbx5^iKO; Nadadur et al., 2016). c-d. Comparison to LA tissue with mild Tbx5 upregulation due to deletion of an intronic regulatory element (Tbx5^Re(int)KO; Bosada et al., 2023). e-f. Comparison to LV tissue with Tbx5 overexpression in cardiomyocytes (Tbx5-OE; this study). g. Comparison of changes in expression of aCM- and vCM-selective genes across all four datasets. Fold-change was calculated between condition with higher Tbx5 (numerator) to condition with lower Tbx5 (denominator). Genes were ordered by ascending Log₂(WT/TBX5^iKO). Genes not detected (ND) for a given experiment are colored white. Normalized enrichment score (NES) and false discovery rate (FDR) values from GSEA using the aCM-selective or vCM-selective gene lists are shown. Negative NES indicates enrichment in the genes upregulated in the lower Tbx5 condition, whereas positive NES indicates enrichment in the genes upregulated in the higher Tbx5 condition.

Extended Data Figure 7.. Characterization of myocyte differentially accessible regions.

Heatmap (a) shows the patterns of differential accessibility between myocyte clusters. The rows contain the union of regions with differential accessibility in the four pairwise comparisons shown in Fig. 5a. ATAC signal in each region is shown for myocyte clusters 1 (control_aCM cluster), 4, 6 (KO_aCM cluster), and 5. The regions are grouped (groups a-f) by their pattern of accessibility change in the four pairwise comparisons. Please also refer to Supp. Table 4. Arrows denote significant enrichment in one cluster compared to another. These six groups fit two predominant patterns: those with and those without predominant accessibility in the control aCM cluster (Myocyte_1). Most regions with predominant accessibility in the control aCM cluster were occupied by TBX5 and had GO terms related to cardiac cell development or striated muscle contraction (b). In contrast, a minority of regions without predominant accessibility in the control aCM cluster were occupied by TBX5 and had GO terms that were atypical for cardiomyocytes (c). Groups with less than 200 regions are not shown in the heatmap.

Extended Data Figure 8.. Proximity and predicted gene linkages demonstrate the regulation of the atrial GRN by control peaks.

a, The different types of association of a genomic region with a gene: (1) Linkage (L). Region-to-gene linkages are predicted based on co-variance of accessibility and expression on a nucleus-by-nucleus basis in the multiome data. (2) Proximity (P). Region-to-gene relationships are inferred by proximity of the region to the gene’s transcriptional start site. (3) Proximity and linkage (P+L). A region-to-gene association can be supported by both proximity and linkage. b, Genes were associated with control and KO regions by linkage (n = 2 control and 2 KO multiome biological replicates). Ratio of gene expression between the control and KO aCM clusters was plotted and compared between groups by the two-sided Mann-Whitney test. ****, P<0.0001. c, Association of control and KO ATAC regions with aCM and vCM genes. Associations were made on the basis of proximity, linkage, or both. d, aCM and vCM genes were interrogated for ass

Extended Data Figure 9.. Motif analysis of control and KO ATAC regions.

a, Aggregation plots for H3K27Ac at control and KO regions. b, Heatmaps of the H3K27Ac signal at control and KO peaks. c, Top-enriched motifs identified in control regions. The TBOX motif was the most enriched, followed by MEF2. An extended table of non-redundant motifs showing the top 21 most enriched motifs in control regions. d, Transcription factor footprinting analysis demonstrates footprints at Tbox and Mef2 motifs in control clusters (Myocyte_1 and Myocyte_4) compared to KO clusters (Myocyte_5 and Myocyte_6). e, Top-enriched motifs identified in KO regions and extended table of the top non-redundant motifs. f. Occupany of control and KO ATAC regions by cardiac TFs in aCMs. Occupancy data is from GEO GSE215065.

Extended Data Figure 10.. Chromatin loops link TBX5 dependent enhancers with atrial genes.

a, Contact maps of Myl7 or Bmp10. Black boxed regions are loops called in each sample and blue boxed regions mark differential loops that are significantly stronger in control samples. b, Genes near control anchors grouped by adjacency to aCM-selective, vCM-selective, or non-chamber selective expression. Control anchors were present near 62 aCM-selective genes and of these, 46 were expressed at greater levels in control Myocyte_1 compared to KO Myocyte_6. Notable genes from previous figures include Nppa, Bmp10, Sbk2 and Myl7. 38 non-chamber selective genes were also upregulated in control samples and linked to enhancers by TBX5-dependent looping. These included Gja1 and Tead1. Only 5 vCM-selective genes were found near control anchors. c, 50 genes neighbored KO anchors. These genes included 5 aCM-selective genes and 3 vCM-selective genes. Most of the differentially expressed genes near KO anchors were more highly expressed in Myocyte_6 (KO) compared to Myocyte_1 (Ctrl).

Figure 1.. Inactivation of Tbx5 in aCMs results in atrial remodeling.

Mice were treated with 2 x10¹¹ VG/g AAV at P2. a, Rosa^mTmG Cre reporter mice were injected with AAV9:Nppa-Cre. Hearts were analyzed at P20. Cre-activated EGFP signal was restricted to the atria. b-d, Whole mount images of control and atrial-specific Tbx5 knockout (Tbx5^AKO) hearts. b-c, brightfield. d, red fluorescence channel. e-g, AAV9:Nppa-Cre mediated recombination. H2B-mCh Cre reporter mice were treated with AAV9:Nppa-Cre. mCherry+ nuclei in atria (e) and ventricles (f) was quantified (g) by confocal imaging of P20 myocardial sections. *, P = 0.012. n=3 hearts per group. h, TBX5 protein level. Atrial and ventricular lysates, prepared at P20, were analyzed by western blotting in biological duplicates. i, RT-qPCR for Tbx5 mRNA from atria and ventricles. **, P = 0.0079; ns, not significant (P=0.48). n=3 control or 4 Tbx5^AKO hearts per group. j-k, Echocardiographic assessment of atrial size. j, representative echo images used to quantify atrial size, compared to aortic size. Bar, 1 mm. j, Relative size of control (n=5) and Tbx5^AKO (n=4) left atria, normalized to aortic diameter. P15, * P = 0.04; P17, * P = 0.011; P20, **** P = 0.0001. Unpaired, two-sided t-tests. Graphs show mean ± SEM. RA, right atrium; LA, left atrium; LV, left ventricle; Ao, aorta.

Figure 2.. Inactivating Tbx5 alters the expression of aCM-selective genes.

a-b, Left atrial myocardial expression of aCM-selective myosin light chain MYL7 (a) or MYL4 (b). Arrows, nuclear mCherry indicative of cell transduction by AAV9:Nppa-Cre. White boxed regions are magnified in insets. Bar: 40 µm. c-d, Expression of vCM-selective myosin light chain MYL2. White boxed regions are enlarged to the right. Arrows, Tbx5^AKO atrial regions that express MYL2. Bar, 100 µm. e, Transcript levels in control and Tbx5^AKO atria. Myl7, Myl4, Myl2, Nppa, and Nppb transcripts were measured from pooled left and right atrial samples of the indicated genotype. Graphs show mean ± SEM. Unpaired two-sided t-test. Myl7: P=2.4E-8, Myl4: P=1.3E-5 Nppa: P=2.5E-6; Myl2, P=0.0041; Nppb, P=0.087. n=9 (Tbx5^AKO) and 7 (control), except for Myl4 control n=8. f, Genome-wide effect of Tbx5 atrial knockout on expression of left atrial genes. Principal component analysis was used to compare transcriptomes of isolated aCMs, isolated vCMs, control LA, and global adult Tbx5 knockout LA (GSE129503). Inactivating Tbx5 shifted aCM gene expression towards vCMs along PC1 (black arrow). g, Definition of aCM- (red) and vCM- (blue) selective transcripts based on P0 RNA-seq of aCMs and vCMs (from GSE215065; Log₂FC > 1.5, P_adj <0.05, Wald test followed by Benjamini and Hochberg correction). h, Enrichment of aCM- and vCM- selective genes defined in (g) among genes down- and up-regulated, respectively, in Tbx5^iKO left atrial samples. Normalized enrichment score (NES) and FDR by Gene Set Enrichment Analysis (GSEA).

Figure 3.. Concurrent scRNA-seq and scATAC-seq analysis of cell states in control and Tbx5^AKO atria.

a, Experimental design. Neonatal Tbx5^Flox/Flox mice were treated with AAV9:Nppa-EGFP (control) or AAV9:Nppa-Cre (KO). At P20, isolated nuclei from two left atria per replicate were analyzed by concurrent scATAC-seq and scRNA-seq. Two replicates were prepared per group. b, Clustering of cell states based on scRNA-seq, scATAC-seq, or both assays using weighted nearest neighbor (WNN) analysis. c, Cluster identities were established using cell-type specific markers. The total number of nuclei and the contribution from control and TBX5^AKO samples are shown. The black dotted line indicates the overall relative contribution of KO and control nuclei to the assay. Other dataset metrics are available in Extended Data Figure 4. CM, cardiomyocyte; Fibro, fibroblast; Peri, pericyte; EC, endothelial cell; Endo, endocardium; Epi, epicardium; Adipo, adipocyte; Macro, macrophage.

Figure 4.. TBX5 is required to promote the expression of aCM genes.

a, Pseudotime trajectory of myocyte clusters. Trajectory analysis was performed with Monocle3. Clusters containing predominantly control (Ctrl) and knockout (AKO) myocytes participated in two separate inferred trajectories. b-c, Biological process GO terms enriched for genes with significantly greater expression in control (pink, Myocyte_1) or KO (magenta, Myocyte_6) clusters. The top 10 terms are shown. Fisher’s Exact test with Benjamini Hochberg multiple testing correction. d, Distribution of chamber-selective genes among genes differentially expressed between control (Myo_1) or KO (Myo_6) myocyte clusters (left column) and from bulk RNA-seq experiment comparing WT and Tbx5^iKO left atria (GSE129503; right column). aCM-selective genes (middle) were disproportionately found in genes with significant upregulation in control aCMs and WT LA, and vCM-selective genes (bottom) were enriched in genes upregulated in KO aCMs or LA. snRNA-seq statistics, Wilcoxon Rank Sum test and Bonferroni correction. Bulk RNA-seq analysis, Wald test and Benjamini Hochberg correction. e, Statistical analysis of distribution of aCM- or vCM-selective genes for the control vs. Tbx5^AKO aCMs (left) and control vs. Tbx5^iKO left atria (right) datasets. Fisher’s exact test (two-tailed).

Figure 5.. Multiomics reveals a set of Tbx5-dependent CREs that promote aCM gene expression.

a, Myocyte trajectories and comparisons performed to identify differentially accessible genomic regions. Arrows indicate the four pairwise comparisons between the clusters at the start and end of the control and Tbx5^AKO trajectories. b, Differentially accessible regions were analyzed for four the pairwise comparisons. The number of differentially accessible regions (region #) is shown at the end of each bar. The bar indicates the fraction of these regions that are occupied by TBX5 in wild-type aCMs. c, Accessibility change between the control and TBX5 KO aCM clusters. Plotting the average log₂FC accessibility of regions in the control aCM cluster compared to the KO aCM cluster revealed 1846 and 309 significant control and KO regions, respectively. (Logistic regression model followed by comparison with a null model using a likelihood ratio test, Bonferonni corrected Padj. < 0.05). Selected genes nearest to top control and KO regions are listed. d, Expression of genes nearest to control or KO regions, from Tbx5^AKO aCM snRNA-seq cluster Myocyte_1 vs Myocyte_6 (left) and Tbx5^iKO LA bulk RNA-seq (right, R26^CreERt2 (n = 4) vs. Tbx5^Flox/Flox;R26^CreERt2 (n = 6)) datasets. Two-sided Mann-Whitney test comparing control and KO ATAC regions. e, GO terms for genes nearest to control regions or KO regions. Genes next to control regions tended to be associated with cardiomyocyte function, while genes adjoining KO regions were associated with TGFβ signaling and remodeling terms (Fisher’s Exact test with Benjamini and Hochberg multiple testing correction).

Figure 6.. TBX5 maintains local chromatin structure to regulate gene expression.

a, Loops identified by H3K27ac HiChIP. Loops were classified as enhancer-enhancer (E-E), enhancer-promoter (E-P), and promoter-promoter (P-P). b, Differential loops between control and Tbx5^AKO aCMs. Statistical analysis identified 510 differential loops with FDR<0.1, 263 with greater loop score in control (“Control loops”) and 247 with greater loop score in Tbx5^AKO (“KO loops”). c, Overlap of differential loop anchors with control and KO ATAC regions. Top row (anchors) considers each loop anchor individually. Bottom row (loops) considers overlap of at least one anchor from each loop. d-e, Comparison of genes associated with control ATAC regions only, control loops only, or both control ATAC regions and control loops. d, Fraction of genes that are aCM-selective. Proportion test vs. “no” association group. e, gene expression ratio in control_aCMs (Myo_1) vs. KO_aCMs (Myo_6). Kruskal-Wallis vs “no” association group. f, Interaction maps for control and Tbx5^AKO samples near the Nppa locus. Significant loops are indicated by boxes on the lower left part of the plot. Interactions significantly enriched in control samples are indicated by blue boxes in the upper right-hand region of the plot. g, Genome browser views of indicated chromatin features. TBX5 bioChIP-seq from aCMs is from GSE215065. Asterisks mark regions with significantly enriched accessibility in control_aCM (Myo_1) compared with KO_aCM (Myo_6).

Figure 7.. TBX5 regulates enhancer activity, accessibility, and looping at aCM-selective genes Myl7 and Bmp10.

a, AAV enhancer-reporter design. Enhancers of Myl7 or Bmp10 (highlighted and labeled in e and f) were cloned within the 3’UTR of the mCherry reporter gene, downstream of a minimal hsp68 promoter. The RNA Pol III U6 promoter drove the expression of a Broccoli non-coding RNA, which was used to normalize for AAV transduction efficiency. b, Activity of Myl7 and Bmp10 enhancers require Tbx5. AAV9 containing either the Myl7 or Bmp10 enhancer was co-injected into P2 Tbx5^Flox/Flox mice with either AAV9:Nppa-EGFP (top) or AAV9:Nppa-Cre (bottom). Hearts were visualized at P8. Number of samples with positive signal and total number of samples is shown. c-d, Quantification of mCherry transcripts. mCherry and Broccoli transcript levels were measured by RT-qPCR from left atrial RNA. mCherry was normalized to Broccoli. Graphs show mean ± SEM. Unpaired two-sided t-test: *, P = 0.0359; ***, P= 0.0005. Bmp10, n = 3 control and 4 Tbx5^AKO. Myl7, n = 4 control and 4 Tbx5^AKO. e-f, Genome browser views of the Myl7 and Bmp10. aCM TBX5 bioChip-seq data is from GSE215065. Statistically significant differences in looping are diagrammed. Asterisks mark regions with significantly reduced accessibility in Tbx5^AKO aCMs.

Figure 8.. TBX5 promotes atrial identity.

a. Identification of aCM- and vCM-selective genes shared between P0 and P28. RNA-seq from P0 and P28 aCMs and vCMs were compared to identify genes that were chamber selective at both stages. b. Enrichment of stage independent aCM genes among genes significantly upregulated in control or KO in the Tbx5^AKO aCMs (left) and Tbx5^iKO left atria (right). c. aCM and vCM identity genes were identified by machine learning. Genes overlapping both cell types were removed, leaving a total of 147 atrial identity (AI) genes and 236 ventricular identity (VI) genes. d. Enrichment of AI and VI genes among genes significantly upregulated in control or KO in the Tbx5^AKO aCMs (left) and Tbx5^iKO left atria (right). e. SEs were identified by H3K27ac signal in aCMs or vCMs. Chamber-selective SEs were enriched for expression of chamber-selective genes. f. Overlap of different classes of SEs with regions with greater accessibility in control (control regions) or Tbx5^AKO (KO regions) aCMs. g. Enrichment of genes neighboring aCM-selective and vCM-selective SEs among genes significantly upregulated in control or KO in the Tbx5^AKO aCMs (left) and Tbx5^iKO left atria (right). b, d, g: Fisher’s exact test, two-tailed.