Single-nucleus multi-omics implicates androgen receptor signaling in cardiomyocytes and NR4A1 regulation in fibroblasts during atrial fibrillation

Abstract

The dysregulation of gene expression programs in the human atria during persistent atrial fibrillation (AF) is not completely understood. Here, we reanalyze bulk RNA-sequencing datasets from two studies (N = 242) and identified 755 differentially expressed genes in left atrial appendages of individuals with persistent AF and non-AF controls. We combined the bulk RNA-sequencing differentially expressed genes with a left atrial appendage single-nucleus multi-omics dataset to assign genes to specific atrial cell types. We found noncoding genes at the IFNG locus (LINC01479, IFNG-AS1) strongly dysregulated in cardiomyocytes. We defined a gene expression signature potentially driven by androgen receptor signaling in cardiomyocytes from individuals with AF. Cell-type-specific gene expression modules suggested an increase in T cell and a decrease in adipocyte and neuronal cell gene expression in AF. Lastly, we showed that reducing NR4A1 expression, a marker of a poorly characterized human atrial fibroblast subtype, fibroblast activation markers, extracellular matrix remodeling and cell proliferation decreased.

Fig. 1. Single-nucleus multi-omics profiling to identify cell types in LAAs.

a, Each cell type is colored differently in the uniform manifold approximation and projection (UMAP). b, The bar plots show the fraction of sample (left) and nucleus count (right) by cell type. Red, AF; Blu, SR. c, Cell-type expression of key marker genes. d, The histogram at the top represents the distribution of ATAC-seq peak count given how many cell types each peak is measured in. The bar plot at the bottom represents the fraction of ATAC-seq peaks located in ENCODE cCREs given how many cell types each peak is measured in. e, Comparison of the ATAC-seq peaks between cell types found in the scAF dataset and the Human Enhancer Atlas. We labeled cell types with the best ATAC-seq peak overlaps. f, Activity of the TF motif (from snATAC-seq, left) and expression of the corresponding TF (from snRNA-seq, right) in each cell type. g, UMAP of motif activity (green) and expression level (red) for two TFs: ESRRG on the left and TCF21 on the right. Adipo, adipocytes; CTCF, CCCTC-binding factor mark; enhD, distal enhancer; enhP, proximal enhancer; K4m3, lysine 4 tri-methyl mark; not.Encode, peak found in the scAF dataset without any overlapping ENCODE cCREs; PC, pericyte; Prom, promoter; SMC, smooth muscle cell.

Fig. 2. A gene expression signature to identify CMs from individuals with AF.

a, To create a CM-specific gene expression signature for AF, we identified genes specifically expressed in CMs in the scAF dataset and selected from those genes that are differentially expressed in two large LAA bulk RNA-seq experiments. This prioritized 62 and 121 genes that are upregulated and downregulated, respectively, in AF. b, For each scAF participant, we plotted the distribution of the upregulated and downregulated expression signatures in CMs. Blue are SR controls and red are individuals with AF. c, We combined CMs from four publicly available cardiac snRNA-seq datasets with scAF. The AF signatures clearly separate CM meta cells from individuals with AF (red) from the rest of the CMs (blue and gray). d, Expression of six CM-specific genes from the AF-upregulated signature in various snRNA-seq datasets. ACM, arrhythmogenic cardiomyopathy; ACM_DCM, arrhythmogenic and dilated cardiomyopathy dataset; Atrial_H.Atlas, atrial heart atlas nuclei dataset; BZ, boarder zone; CTRL_A/D, control samples from the ACM_DCM dataset; CTRL_MI, control samples from the MI dataset; DCM, dilated cardiomyopathy; FZ, fibrotic zone; IZ, ischemic zone; LA, left atria; LAA_AF, this study scAF dataset; MI, myocardial infarction dataset; NCM, non-compaction cardiomyopathy; RA, right atria; RZ, remote zone. e, The plot at the top shows the correlation coefficients (Pearson’s r) between the AF-upregulated signature and TF motif activities in CMs. The scatterplot at the bottom captures the correlation between TF motif and the AF-upregulated signature (x axis) and between TF expression level and the AF-upregulated signature in CMs from the scAF dataset (y axis). Red points highlight TFs with significant correlations at the level of the motif and expression with the AF-upregulated signature (FDR < 0.01). f, In the scAF CMs, the AF-upregulated signature is correlated with the AR motif activity (snRNA-seq; top) and with AR expression (middle). The AR motif and expression are also correlated in CMs (bottom). Blue indicates SR controls and red indicates individuals with AF. For each plot, we calculated Pearson’s r and a nominal two-sided P value. g, The Tn5 footprinting shows a stronger enrichment near AR motifs in CMs versus other cell types, and in CMs from individuals with AF versus SR controls. The results by cell type are shown in Supplementary Fig. 11.

Fig. 3. FB subtypes in human LAAs.

a–d, Different analyses of FBs in our scAF multiome dataset. e–g, Summary of analyses in which the scAF FB was integrated with FBs from additional cardiac single-nucleus datasets. a, FB UMAP colored by FB subtype. b, Dot plot showing averaged normalized expressions for the strongest marker genes for each FB subtype. c, Bar plot showing the FB subtype proportion by rhythm. d, Dot plots showing the top three gene sets (by enrichment scores) in each FB subtype from a gene-set overrepresentation analysis based on specificity of expression. For this analysis, we used the Hallmark and Gene Ontology Biological Process (GO BP) gene-set libraries. Gene-set P values were adjusted (P_adj) for multiple testing (Benjamini–Hochberg method) below 0.05 (fast gene-set enrichment analysis enrichment P values are estimated using an adaptive multilevel split Monte Carlo method). e,f, Integrated FB UMAP of four human cardiac datasets (Methods). e, Eight clusters labeled by their marker genes (the ‘Not.specific’ cluster designates a small proportion of cells for which we could not identify a strong gene marker). f, UMAP positions of the aFB3 cluster from the scAF dataset. g, Corresponding fractions of the FB subtypes identified in scAF to FB subtypes identified in the multi-study integrated FB clusters in e. NES, normalized enrichment score; NF-κB, nuclear factor kappa B.

Fig. 4. NR4A1 deficiency reduces expression of fibrotic markers and alters the function of human atrial FBs.

a, NR4A1 mRNA levels (quantified by quantitative PCR with reverse transcription (RT–qPCR), relative to the geometric mean of housekeeping genes HPRT1 and RS18) in human atrial FBs. b, Representative immunoblots of NR4A1 protein (normalized to GAPDH) in human atrial FBs. c,d, siRNA-mediated knockdown of NR4A1 (NR4A1-KD) mRNA (analyzed by RT–qPCR) and protein (immunoblotting) in human atrial FBs. e–n, Effect of NR4A1-KD on expression of selected fibrotic markers: αSMA (e and f), periostin (g and h), fibronectin (i and j), type I collagen (Col 1; k and l) and type III collagen (m and n) in atrial FBs. o, Effects of NR4A1-KD on 24-h proliferation of human atrial FBs. p, Effects of NR4A1-KD on cell migration (scratch assay) at 8-h and 12-h time points. Scale bar, 200 µm. Data are the mean ± s.e.m. (a, b and p) or mean with paired lines (c–o); each dot (a and b) or paired line (c–o) represents an independent biological replicate. Data are expressed as a percentage of the control group (NC-siRNA; c–o) or the percentage at 0 h (p). P values were determined from raw values (2^−∆Ct) and reported as two sided using a paired t-test (c, d, j and o), ratio paired t-test (e–i, k, m and n), Wilcoxon test (l) and two-way analysis of variance (ANOVA) with Sidak correction (p). N, individual biological donors. Source data

Extended Data Fig. 1. Characterization of two left atrial appendage (LAA) bulk RNAseq dataset.

The left and right panels show results from the CTSN cohort (31 AF and 31 SR samples) and GSE69890 (130 AF and 50 SR samples), respectively. (a) Principal component analysis (PCA) of the 500 most variable genes in each dataset. (b) Scatter plot showing the log₂-transformed expression of the epicardium marker ITLN1 against the first principal component (PC1). (c) Top ten genes with highest loadings for PC1. (d) Dotplot showing average normalized expressions for genes in (c) by cell-type in the scAF single-nucleus RNAseq LAA dataset. Adipo; Adipocytes, CM; Cardiomyocytes, EC; Endothelial cells, FB; Fibroblasts, PC; Pericytes, SMC; Smooth muscle cells.

Extended Data Fig. 2. Impact of stand information of differential gene expression in bulk left atrial appendage (LAA) RNAseq datasets.

(a) Miami plot comparing the bulk RNAseq differential expression between atrial fibrillation (AF) and sinus rhythm (SR) samples -log10(false discovery rates [FDR]) in the CTSN and GSE69890 datasets. Colors and numbers on the x-axis show the chromosomal position of each gene. Red crossed genes denote genes that were only differentially expressed when omitting strand information in the CTSN dataset (Supplementary Table 3). The dashed red box highlights the IFNG locus. (b) Gene expression profiles at the IFNG locus in the CTSN cohort. (top) Coverage plot of all AF and SR samples combined. (center) -log10(FDR) with (orange) and without (blue) strand information provided during read pseudoalignment (Methods). Signs in parentheses show the strand information for each gene. (bottom) Transcript annotation at this locus colored by strand.

Extended Data Fig. 3. Integration of bulk and single-nucleus RNAseq data in left atrial appendages (LAAs) to identify differentially expressed genes (DEGs) between patients with atrial fibrillation (AF) or control in sinus rhythm (SR).

(a) Signed log₁₀(p-value) for pseudo-bulk differential gene expression analysis of each cell-type between AF and SR samples. The sign is based on the direction of the log₂ fold-change. Gene colored in red are differentially expressed (false discovery rate [FDR] <0.05). We labeled genes that were also differentially expressed with concordant direction of effect in the two bulk RNAseq datasets. The inset histograms show the number of (top) upregulated and (bottom) downregulated DEGs in each cell-type. (b) Volcano plot showing the results of differential expression analysis in the CTSN cohort (LAA bulk RNAseq dataset). Genes colored in red have an FDR < 0.05 and an |log₂ fold-change | >0.25. (c) Single-nucleus uniform manifold approximation and projection (UMAP) colored based on the level of normalized expression of (top) LINC01479 and (bottom) IFNG-AS1 in (left) AF and (right) SR patients. This shows increased expression of both genes (blue) in cardiomyocytes of AF patients. (d) Violin plots showing the distribution of the cardiomyocytes normalized expression of (top) LINC01479 and (bottom) IFNG-AS1 in each LAA sample analyzed in the scAF dataset. Violin distributions of samples colored in red and blue are AF and SR, respectively. (e) Boxplots showing the normalized counts of (left) LINC01479 and (right) IFNG-AS1 in the two bulk RNAseq cohorts used in this study. J. Hsu corresponds to the NCBI GEO dataset GSE69890. The box plots show the upper quartile, median and lower quartile. Adipo; Adipocytes, CM; Cardiomyocytes, EC; Endothelial cells, FB; Fibroblasts, PC; Pericytes, SMC; Smooth muscle cells.

Extended Data Fig. 4. Gene expression module analysis using WGCNA.

(a) Dendrogram of gene x gene similarity based on their expression in left atrial appendage bulk RNAseq datasets. Each gene’s module attribution by WGCNA is labeled by colors at the bottom. (b) Heatmap shows the scaled mean module score (centered and divided by the module standard deviation) by cell-type in the scAF dataset. The dendrograms show similarity in terms of module (x-axis) and cell-type (y-axis). Adipo; Adipocytes, CM; Cardiomyocytes, EC; Endothelial cells, FB; Fibroblasts, PC; Pericytes, SMC; Smooth muscle cells.

Extended Data Fig. 5. WCGNA bulk RNAseq module gene set overrepresentation analysis.

(left) Dot plots that show the top three gene sets from an overrepresentation analysis of the differentially expressed genes (DEGs) in each module. For this analysis, we used the PanglaoDB and gene ontology biological process (GO BP) gene set libraries. (center) Volcano plots that shows the log₂ fold-change and -log₁₀(false discovery rate [FDR]) statistics from the DEG analysis, stratified by module. Red and blue integers indicate the number of atrial fibrillation upregulated and downregulated genes, respectively, in each module. (right) Violin plots showing the distribution of the module scores in each cell-type of the scAF dataset. Adipo; Adipocytes, CM; Cardiomyocytes, EC; Endothelial cells, FB; Fibroblasts, PC; Pericytes, SMC; Smooth muscle cells.

Extended Data Fig. 6. Subcluster analysis of cardiomyocytes (CMs) in the scAF single-nucleus RNAseq dataset did not identify atrial fibrillation (AF)-specific cell clusters, but rather an artifact due to a pre-existing additional heart condition in one of the donors.

(a) CM uniform manifold approximation and projection (UMAP) colored by the two clusters identified by Seurat. (b) Dotplot showing average normalized expressions for the top six marker genes for each CM subcluster in a. (c) Barplot showing the sample portions in each CM subcluster. Red and blue samples represent AF and sinus rhythm (SR) patients, respectively. (d) CM UMAP as in a with each donor sample colored differently. c and d indicate that most of the CMs in subcluster 1 belong to donor CF91. This participant was the only one who suffered from a myocardial infarction less than a year before tissue collection, possibly explaining the strong difference of its CM transcriptome. (e-f) Scatter plot of (e) NR3C2 motif activity and (f) expression correlations with the AF signature UP scores in CM metacells. Red and blue samples represent AF and SR patients respectively. (g) Ventricular CM gene scores derived from 3 CM states from normal and myocardial infarction human hearts (Methods) in CM from the scAF dataset.

Extended Data Fig. 7. Testing the AF UP and DOWN gene expression signatures in cardiomyocytes (CMs) from independent single-nucleus RNAseq datasets.

(a) Barplot showing the number of samples found in each dataset used to compare the specificity of the AF signatures. (b) Violin plots showing the distribution of the AF signature scores of CM metacells in each disease or cardiac chamber. All comparisons with the scAF AF distribution (leftmost) are significant (t-test P < 2.2×10^-16, with the caveat that meta-cells from the same donor are not independent). ACM; arrhythmogenic cardiomyopathy, ACM_DCM; arrhythmogenic and dilated cardiomyopathy dataset, AF; atrial fibrillation, Atrial_H.Atlas; atrial heart atlas nuclei dataset, BZ; boarder zone, CTRL_A/D; control samples from the ACM_DCM dataset, CTRL_MI; control samples from the MI dataset, DCM; dilated cardiomyopathy, FZ; fibrotic zone, IZ; ischemic zone, LA; left atria, LAA_AF; this study scAF dataset, MI; myocardial infarction dataset, NCM; non-compaction cardiomyopathy, RA; right atria, RZ; remote zone, SR; sinus rhythm.

Extended Data Fig. 8. Comparison of the atrial fibrillation (AF) UP and DOWN gene expression signatures when using our scAF and the single-nucleus RNAseq data from the Heart Atlas to select cardiomyocyte (CM)-specific genes.

To create the signatures, we intersected genes that are differentially expressed between atrial fibrillation (AF) and sinus rhythm (SR) left atrial appendages in two bulk RNAseq experiments, together with genes that are specifically expressed in CMs as determined by single-nucleus RNAseq using specificity metrics (area under the curve >0.5, Methods). The results presented here show that very concordant AF signatures are obtained when using either our scAF data or the single-nucleus RNAseq data from the Heart Atlas, thus alleviating concerns of data overfitting.

Extended Data Fig. 9. Validation of aFB3 in external single-nucleus RNAseq of left atrial appendages obtained from patients with AF or in normal SR.

(a-b) Uniform manifold approximation and projection (UMAP) of the fibroblasts (FB) in the GSE224959 single-nucleus RNAseq dataset, showing (a) FB subtype and (b) NR4A1 density of expression. (c) Scatter plot of gene specificity of expression measured by area under the receiving operating curve (AUC) for aFB3 of the scAF dataset and cluster 3 in (a). The Pearson’s p-value and correlation coefficient are shown.

Extended Data Fig. 10. Sub-analysis showed sex-specific differences in NR4A1 regulation of fibrotic responses in human atrial fibroblasts (FBs).

(a-k) Effect of NR4A1 siRNA-mediated knockdown (NR4A1-KD) on expression of selected fibrotic markers: αSMA (a-b), periostin (c-d), fibronectin (e-f), collagen 1 (g-h) and collagen 3 (i-j) in male atrial FBs (blue bars). (k) Effects of NR4A1-KD on 24-hour proliferation of male FBs. (l-v) Effect of NR4A1-KD on expression of selected fibrotic markers: αSMA (l-m), periostin (n-o), fibronectin (p-q), collagen 1 (r-s) and collagen 3 (t-u) in female atrial FBs (red bars). (v) Effects of NR4A1-KD on 24-hour proliferation of female FBs. Data are mean with paired lines; each paired line represents an independent biological replicate. Expression is expressed as % of control group (NC-siRNA). P values were determined from raw values (2^-∆Ct) using ratio paired t-test (a-g, i-j, l-u) and paired t-test (h, k, v). N, individual biological donors. Source data