Dissecting regulatory non-coding GWAS loci reveals fibroblast causal genes with pathophysiological relevance to heart failure

Abstract

Heart failure is caused in part by cardiac remodeling processes that include the death of cardiac myocytes and their replacement by cardiac fibroblasts. Here, we hypothesize that cardiac fibroblasts may harbor epigenetic contexts in which heart disease-associated non-coding SNPs perturb gene expression relevant to disease. To test this, we utilized male primary cardiac fibroblasts to generate high-resolution Hi-C data and integrate it with functional genomic information to annotate and link putative distal regulatory elements in heart disease-associated loci to gene promoters. We identify several target genes with established roles in cardiac fibrosis and/or heart disease (GJA1, TBC1D32, CXCL12, IL6R, and FURIN). We perform Perturb-seq in immortalized male cardiac fibroblasts to knock out putative regulatory elements, confirming regulatory relationships involving GJA1, CXCL12, and FURIN. Our results demonstrate that multi-omic approaches can delineate pathophysiologically relevant regulatory circuits connecting protein-coding genes to non-coding genetic variants associated with disease.

Heart failure can be caused by cardiac fibroblasts replacing myocytes. Here, the authors use functional genomic data from fibroblasts, genetic signals enriched in people with heart disease, and gene perturbation analyses to link disease-associated regulatory elements to protein-coding genes.

Introduction

Heart failure is characterized by systolic and/or diastolic dysfunction and represents a rapidly growing public health problem and unmet medical need in the US. While the population prevalence of heart failure was 2.4% in 2012, this figure will increase more than ten-fold to 25% by 2030 as a surplus of Americans reach age 65¹. In parallel, the medical costs of heart failure are projected to rise from $20.9 billion in 2012 to $53.1 billion in 2030². Heart failure is a complex syndrome secondary to functional or structural deficits that impair ventricular filling or blood ejection, including ischemic and valvular heart disease, cardiomyopathy, congenital heart disease, and pericardial disease. The diminished contractility and reduced cardiac output due to these initiating insults and cardiac myocyte death lead to maladaptive compensatory cardiac remodeling (hypertrophy or dilation concurrent with cardiac fibrosis); this stiffens the myocardium and further impairs cardiac filling and myocardial contraction³, ultimately leading to cardiac decompensation and heart failure.

Cardiac fibroblasts are highly relevant to the pathophysiology of heart failure and its pathogenesis. These mesenchymal cells are involved in remodeling-associated hypertrophy and dilation, as well as post-infarction replacement of dead cardiomyocytes with a collagen-rich (fibrotic) scar to maintain the heart’s structural integrity. However, if disease-activated fibroblasts become hyperactive (myofibroblasts), their excessive production and secretion of extracellular matrix (ECM) can result in the appearance of patchy or diffuse excess fibrotic tissue, and, consequently, increased myocardial wall stiffness. Fibroblast proliferation and phenoconversion to myofibroblasts are both initiated by TGF-β and its protease regulators such as FURIN⁴ and MMP family proteins⁵. This profibrotic signal is mediated by transcription factors such as SMAD3 and SMAD4, as well as CCN2/CTGF⁶, and is further modulated by KLF15⁷. Other important factors in cardiac fibrosis include Ang II⁸, EDN1^{9, 10–11}, Wnt-β-Catenin signaling^9,12, and the SDF1α (CXCL12)/CXCR4/CXCR7 axis¹³.

Heart failure and its antecedent diseases have a strong genetic component and have accordingly been intensively investigated in genome-wide association studies (GWAS), which have yielded numerous loci with robust statistical associations. However, across cardiovascular and many other unrelated diseases, most associated SNPs reside within non-coding (intronic or intergenic) parts of the genome. While most ( ~80%) of these loci are themselves enriched in putative regulatory elements^14,15, their target genes are frequently not known. One approach to assign target genes to GWAS variants is expression quantitative trait locus (eQTL) mapping. However, the limited utility of eQTL overlap alone for understanding GWAS is limited: most existing eQTL signals have been generated from whole tissues comprised of heterogeneous mixtures of cells and discovered in small cohorts, and consequently subject to different linkage disequilibrium patterns and potential confounders^16,17. Therefore, even when co-localized with eQTLs, non-coding GWAS signals are often mapped, and mis-attributed, to the nearest gene(s), and their direct relevance to disease is not clear.

Studies in numerous disease contexts suggest that causal SNPs’ regulatory targets may often be thousands to millions of base pairs away, yet physically proximal due to individual DNA loops^{18, 19, 20, 21–22}, as well as larger-scale topologically associating domains (TADs)^{23, 24–25}. Recent advances in chromosome conformation capture techniques (i.e. Hi-C) now permit genome-wide identification of significant DNA-DNA interactions and TADs. However, the majority of published Hi-C data are low-resolution²⁶ or derived from cell lines not representative of disease contexts²⁷. Capture Hi-C partially allays resolution issues by identifying loops involving promoter-centric regions at high-resolution²⁸, but in doing so may introduce capture efficiency biases²⁹.

Previous studies have applied capture Hi-C or low-resolution (5-kb) Hi-C in cardiomyocytes derived from stem cells or whole tissue^{30, 31, 32–33}; in addition to problems associated with biased interaction sampling, the data generated in these studies had little to no cardiac fibroblast representation. Notably a mouse study that assayed genome-wide single cell open chromatin in whole hearts found enrichment of heart failure-associated genetic signals in fibroblasts³⁴, further affirming the importance of this cell type to heart failure.

To address these gaps and determine how non-coding GWAS loci contribute to heart failure etiology, we integrated high-resolution 3D chromatin interaction and open chromatin data in primary human ventricular fibroblasts. Using this multi-layered omics dataset, we systematically constructed genome-wide gene regulatory circuits between heart disease GWAS SNPs and their target genes, irrespective of genetic distance. Significant Hi-C loops were intersected with pairs of known gene promoters and heart disease loci prioritized with ATAC-seq and RNA-seq data to identify regulatory relationships and contrasted with GTEx eQTLs. The resultant regulatory circuits’ target genes were enriched for established cardiac remodeling and heart disease associated genes and may include potential targets for therapeutic intervention. Finally, we validated a subset of high priority SNP-gene pairs using CRISPR-Cas9 deletion of the regulatory regions followed by single cell RNA-seq, confirming their regulatory role and the value of multi-omic data integration for target gene identification and prioritization.

Results

Characterization of functional genomic features relevant to heart failure and gene regulation

Given the critical roles of cardiac fibroblasts in the pathogenesis of heart failure, we hypothesized that the integration of chromatin interactions from these cells with heart disease-associated loci would reveal gene regulatory circuitry underlying heart failure. To capture chromatin interactions, we deeply sequenced large Hi-C libraries to identify long- and short-range interactions across the genome in primary Normal Human Cardiac Fibroblasts – Ventricular (NHCFV) cells. NHCFV cells are commercially available primary cells derived from non-diseased human tissues.

The HiC-Pro and Juicer pipelines both called nearly 1.5 billion valid interaction pairs after removing duplicate pairs ( ~20%). Over 1 billion (~85%) of these interactions were in cis (on the same chromosome) and a large proportion (60–66%) were long-range ( ≥20-kb) (Table S1), indicative of high-quality library construction²⁷. Our Hi-C data had 2-kb resolution, which should be sufficient to unambiguously capture interactions involving most enhancers (mean length 540.7 ± 289.8 bp, and up to 1792 bp) and super-enhancers (mean length 3324.5 ± 833.7 bp, with a maximum size of 6345 bp), based on FANTOM5 annotations³⁵. Having identified the scaffolding underlying gene regulation in NHCFV cells, we proceeded to annotate the cell type-specific pairwise interactions with functional elements.

A subset of the numerous open chromatin regions within cells represents sites of transcriptional regulation during homeostatic conditions, and their perturbation by disease-associated genetic variation may modulate gene expression during disease. To define these putative regulatory elements in cardiac fibroblasts, we performed ATAC-seq in NHCFV and immortalized human cardiac fibroblast (iHCF) cells (Table S2). Further analysis of these euchromatic regions identified predicted transcription factor binding sites, which provided additional evidence that these elements participate in gene regulation. The ENCODE ATAC-seq pipeline identified over 250,000 peaks that were shared by pseudo-replicates (naïve overlap peaks) derived from single NHCFV and iHCF ATAC-seq libraries.

HMM-based IdeNtification of Transcription factor footprints (HINT) identified over 300,000 footprinted regions each in NHCFV and iHCF cells. After merging overlapping footprinted regions (footprints containing multiple predicted TF binding sites) within each cell type, we found that 67.6% of unique footprinted regions were common between the primary and immortalized cells. Of these shared regions, 157,307 (nearly 40%) also spanned open chromatin regions. The footprinted regions contained 2,694,471 unique motifs corresponding to binding sites corresponding to 840 known TFs in HOCOMOCO and JASPAR. Notably, several well-characterized cardiac-relevant transcription factor binding sites, including those of TEAD and SMAD family proteins, were enriched in open chromatin regions involved in interactions with accessible promoters in NHCFV cells (Supplementary Data 1).

We also performed bulk RNA-seq on NHCFV, iHCF, and primary human cardiomyocyte (HCM) cells to further characterize these cells, as well as to prioritize chromatin interactions by quantifying gene expression in our cell types of interest (Figure S1 and Table S3). The expression distribution (log₂FPKQ) of expressed genes and Pearson correlation coefficient (r) reveal that iHCF are more similar to NHCFV than HCMs, as expected (Figure S1A). NHCFV and iHCF shared expression of cardiac fibroblast genes (e.g. ACTA2, AIFM2), and low expression for cardiomyocyte markers (e.g. TTN, MEF2C), though we observed differences in certain fibroblast marker genes (e.g. DCN), which is expected when comparing primary and immortalized cells (Figure S1B). Overall, this transcriptomic data layer informs which genes identified from primary fibroblasts can be tested in a cell line that is amenable to genetic perturbation.

Data integration to generate regulatory circuits

We integrated multiple layers of functional genomic data to identify interacting pairs of SNPs (perturbations in potential disease-related regulatory regions) and TSSs (target genes, whose dysregulation may underlie GWAS signals). Figures 1A, B summarize our analysis workflow. We intersected significantly interacting bin pairs from NHCFV cells (Fig. 1C) with accessible SNPs and TSSs to obtain accessible SNP-TSS pairs, and further identified pairs with SNPs in footprints of expressed TFs as well as TSSs of expressed genes. We then annotated SNPs that were in linkage disequilibrium (LD) with GTEx tissue eQTLs associated with genes corresponding to their interacting TSSs.

Fig. 1 Tracing regulatory circuits from genome-wide association study (GWAS) signals to target genes using functional genomics. [Images not available. See PDF.]

A Graphical summary of analyses annotating GWAS loci with 3D chromatin interaction and open chromatin data to prioritize functional variants and identify target genes that may underlie disease associations (P < 5E-8 for one-sided Chi-square tests of association between single nucleotide polymorphism (SNP) genotypes and phenotypes). Significant eQTLs (FDR < 0.05) further corroborated these SNP-gene relationships. B Analysis workflow for integration of genetic, expression, and functional genomic data. We layered significant pairwise interactions (chromatin loops) with SNPs and transcription start sites (TSSs) filtered for assay for transposase-accessible chromatin (ATAC)-seq regions and highlighted those in expressed transcription factors (TFs) footprints (from the HOCOMOCO and JASPAR databases) and expressed target genes, as well as Genotype Tissue Expression Project (GTEx) expression quantitative trait loci (eQTLs), to identify putative active and functional SNP-gene pairs. C Genome-wide contact maps of high-resolution Hi-C interactions from primary normal human cardiac fibroblasts –ventricular (NHCFV). (D) Number of GWAS SNPs associated with heart disease-related traits, including heart failure.

Our comprehensive search for genetic variants relevant to heart failure yielded 5534 SNP-phenotype associations, representing 3880 unique SNPs across heart disease-related pathologies (Fig. 1D). Our search utilized broad criteria to include heart failure risk factors (e.g. atrial fibrillation, coronary artery disease, hypertension) as well as diseases related to metabolic syndrome. This approach enabled identification of loci linked to multiple cardiometabolic traits, both related and unrelated to cardiac remodeling and fibrosis in heart failure. We augmented the heart disease SNPs with proxy SNPs (r² > 0.8 in the European ancestry population), resulting in a total of 73,203 variants, representing 3090 interconnected clusters of correlated SNP pairs (“loci”).

We next proceeded to integrate heart disease-associated genetic variants and gene promoters with ATAC-seq, Hi-C, and RNA-seq data in NHCFV cells (Fig. 2A). First, we intersected open chromatin regions derived from ATAC-seq in NHCFV cells with heart disease-related SNPs and genome-wide TSSs. About 15% of accessible SNPs overlapped with HINT footprint regions, and a smaller subset (7-9%) overlapped with TFBSs of known TFs from HOCOMOCO and JASPAR (Fig. 2B, and Table S4). Most accessible SNPs and TSSs overlapped with significant interactions (chromatin loops) called from our Hi-C data, although the degree of this overlap varied between GOTHiC and HiCCUPS (Table S5). These interactions were all intra-chromosomal and predominantly long-range ( ≥ 20-kb). While GOTHiC called more short-range interactions overall, HiCCUPs primarily identified long-range SNP-gene pairs (Table S6A). We used ATAC-seq and bulk RNA-seq to identify loops involving SNPs in putative active enhancers (accessible to transcriptional machinery) and active genes’ promoters (accessible TSSs of expressed genes) (their properties are described in Table S6B). Although putative enhancers containing footprints with expressed TFs were more likely to be functional, we did not hard-filter for this attribute, in order to be inclusive of regions that may contain TF binding motifs that have not been discovered. Using bulk RNA-seq, several candidate interactions were excluded despite their overlap with putative enhancers overlapping GWAS loci, due to almost undetectable expression of predicted target genes (e.g. GLP2R, NAALAD, CPNE5). In total, we identified over 3000 unique SNP-gene pairs in each cell type using GOTHiC and HiCCUPS. There was a modest overlap between SNP-gene pairs called by GOTHiC (short-range) and HiCCUPS (long-range) (Figure S2), and GOTHiC identified the same SNP-gene pairs in both cell types more frequently compared to HiCCUPS. Notably, over two-thirds of SNP-gene pairs involved target genes that were distal to the SNP reported in the original GWAS, and therefore involved target genes that were not previously reported in GWAS. This proportion is also consistent with the overabundance of long-range cis interactions. A small subset of these SNPs was in LD with significant GTEx eQTLs associated with their target genes in heart regions ( ~2%), although a larger proportion ( ~4%) corresponded to significant GTEx eQTLs in any human tissue (Supplemental Notes).

Fig. 2 Prioritization of active and functional heart disease-associated SNP-gene pairs in normal human cardiac fibroblasts – ventricular (NHCFV) combining genetics, epigenetics, and transcriptomics. [Images not available. See PDF.]

A Sankey diagram of the cell type-specific variation observed in heart disease-associated SNPs and gene accessibility after filtering for open chromatin, chromatin interaction, and expression features, as outlined in Fig. 1B. B Quantification of accessible SNPs identified in NHCFV that are within known transcription factor (TF) footprints of TFs that also have detectable mRNA expression. TADs topologically associating domains; ATAC assay for transposase-accessible chromatin; SNP single nucleotide polymorphism; LD linkage disequilibrium; CEU Utah residents (Center d′Etudes du Polymorphisme Humain) with Northern and Western European ancestry; GWAS genome-wide association study.

SNP-Gene Interactions Involving Known Heart Disease Genes

For a subset of the significant SNP-gene interactions, our multi-omic approach revealed high concordance between chromatin regulation, gene expression, and involvement of TFs and target genes with established roles in the pathophysiology of heart disease, cardiac fibrosis, and heart failure (Table 1). We deeply analytically and experimentally characterize these four gene regulatory relationships and also provide the full list of prioritized putative enhancer-promoter pairs in Supplementary Data 2. These gene regulatory networks reveal epigenetic contexts that underlie human genetics-driven disease-modifying changes in cellular behavior.

Table 1. Significant interactions between heart disease GWAS SNPs involving known heart disease genes and/or transcription factors

SNP	SNP Position	Associated Diseases	Mapped Genes	Min GWAS p value	Min GTEx eQTL p value	TF Footprints (motif match)	Gene	Gene Expression (FPKQ) in NHCFV, HCM, and iHCF cells	TSS(s)	TSS Range	Loop Source	Loop Bin Distance	Loop Reads	Min Loop q-value
rs1919875	chr6: 121,849,906	AF	AL603865.1, HMGB3P18, HSF2	3.0E-11	N/A	N/A	GJA1	325.5 40.7	P2	chr6:121,435,591-121,435,766	GOTHiC	414,000	14	3.28E-22
rs7772537	chr6: 122,185,395	AF	GJA1, LOC105377979, RNU1-18P	1.0E-08	N/A	FOXJ3 (13.7)	TBC1D32	1.0 0.1	P1, P2, P3	chr6:121,334,406-121,334,739	GOTHiC HiCCUPS	850,000	6-68	9.15E-19
rs1680636	chr10:44,180,492	CAD	CXCL12, LINC00841	2.0E-16	N/A	ARNT (11.3) ARNTL (14.1) ATF3 (10.8) BHLHE40 (13.0) BHLHE41 (11.6) CREB3L2 (12.1) MLX (13.7) SREBF1 (11.0) SREBF2 (9.8) TFE3 (14.1) TFEB (12.6)	CXCL12	145.0 498.4	P1	chr10:44,376,910-44,385,105	HiCCUPS	200,000	58	5.84E-04
rs7549250	chr1:154,431,860	MI, CAD	IL6R	6.0E-16	1.6E-20	KLF12 (10.2)	IL6R	0.4 0.3	P19	chr1:154,405,199-154,434,982	GOTHiC	24,000	11	2.14E-84-1.18E-24
rs4129267	chr1:154,453,787	C-reactive protein		2.0E-14	2.4E-09	ELK3 (11.1) REL (13.6)			P18			18,000	21	1.61E-54
rs35346340	chr15:90,884,641	MI, CAD	FURIN, FES	2.0E-25	1.1E-15	EGR1 (11.9) EGR2 (11.9) KLF15 (11.3) KLF6 (10.5) MAZ (13.1) MESP1 (12.8) MXI1 (10.9) SNAI1 (12.0) SP1 (9.9) SP2 (10.1) SP3 (12.2) SP4 (10.2) VEZF1 (10.2) ZEB1 (11.7) ZFX (11.1)	FURIN	77.7 65.1	P4, P6, P10, P11	chr15:90,868,580-90,875,744	GOTHiC	12,000	12	1.21E-24

Transcription factors involved in cardiac fibrosis and/or heart disease are in bold.

SNP single nucleotide polymorphism, GWAS genome-wide association study, GTEx Genotype Tissue Expression Project, TF transcription factor, FPKQ fragments per kilobase per million sequenced, quantile normalized; TSS transcription start site

AF Atrial fibrillation; CAD coronary artery disease; MI myocardial infarction.

NHCFV Normal human cardiac fibroblasts (ventricular); iHCF Immortalized human cardiac fibroblasts.

GWAS p values were obtained from one-sided Chi-square tests of association between SNP genotypes and phenotypes, with genome-wide multiple testing-corrected alpha values of 5E-8.

eQTL q values were obtained from linear regression of gene expression on SNP genotypes, and significant associations were identified based on FDR q values < 0.05.

Significant Hi-C interactions (loops) between binned regions of the genome were based on FDR < 0.1 calculated by the HiCCUPS algorithm. The HiCCUPS statistical test corresponds to a one-sided Poisson test comparing observed vs. expected contacts in four local neighborhoods (donut, horizontal, vertical, and lower-left). BH-FDR multiple testing correction is applied within lambda-bins. In addition to HiCCUPS, we used GOTHiC to identify Hi-C loops. GOTHiC loop testing is performed via a one-sided cumulative binomial test for enrichment, which is again corrected for multiple testing using BH-FDR.

rs7772537 is an intergenic SNP that is in LD with rs12664873, which was identified in a GWAS of atrial fibrillation (OR 1.08 [1.05-1.12], P = 1.2E-8)³⁶ (Table 1, Fig. 3). rs7772537 resides within a TF footprint of FOXJ3 in cardiac fibroblasts and physically interacts with the promoter of TBC1D32 over 850-kb away. An additional independent atrial fibrillation signal was reported at rs1919875, a variant between rs7772537 and TBC1D32 (OR 1.08, P = 3.0E-11)³⁷. rs1919875 is in an open chromatin region and tags the intergenic variant rs13195459 (mapped to HSF2) but interacts with GJA1 over 400-kb away. The rs7772537-TBC1D32 and rs1919875-GJA1 loops are proximal to TAD boundaries called by all four algorithms, particularly Arrowhead and Insulation Score (Fig. 3, and Figure S3). Notably, rs7772537 and rs1919875 do not tag GTEx eQTLs for TBC1D32 or GJA1. GJA1 is highly expressed in NHCFV and iHCF cells, while TBC1D32 is lowly expressed in these cells.

Fig. 3 Significant interactions between an atrial fibrillation-associated putative regulatory elements, at rs7772537 and rs1919875, with TBC1D32 and GJA1, respectively, in cardiac fibroblasts, overlaid with GRiNCH TAD calls. [Images not available. See PDF.]

Diamonds in Hi-C contact maps represent significant interactions (FDR < 0.1) between binned regions of the genome based on the HiCCUPS algorithm. The HiCCUPS statistical test corresponds to a one-sided Poisson test comparing observed vs. expected contacts in four local neighborhoods (donut, horizontal, vertical, and lower-left). BH-FDR multiple testing correction is applied within lambda-bins. Dotted lines in Manhattan plots indicate genome-wide statistical significance thresholds of 5E-8 for one-sided Chi-square tests of association between SNP genotypes and phenotypes. TADs topologically associating domains, ATAC assay for transposase-accessible chromatin, SNP single nucleotide polymorphism, LD linkage disequilibrium; CEU Utah residents (Center d′Etudes du Polymorphisme Humain) with Northern and Western European ancestry; GWAS genome-wide association study.

rs1680636 is an intergenic SNP in footprints of several TFs and interacts exclusively with CXCL12 in a 200-kb loop, in cardiac fibroblasts (Table 1, Fig. 4). CXCL12 is highly expressed in NHCFV and iHCF cells. rs1680636 is an intergenic SNP mapping to both CXCL12 and LINC00841, and is a tag SNP for the intergenic variant rs2457480 associated with CAD (OR = 1.111 [1.109-1.115], P = 2.0E-16)³⁸ (Table 1, Fig. 4). Consistent with the abundance of interactions at the locus surrounding CXCL12, boundaries of larger TADs were identified across most TAD calling tools (Fig. 4, and Figure S4). CXCL12, whose serum levels have been linked to atherosclerosis and atrial fibrillation³⁹, is highly expressed in NHCFV and iHCF cells, but rs1680636 and its proxies are not CXCL12 eQTLs in GTEx tissues.

Fig. 4 Significant interaction between a coronary artery disease-associated open chromatin at rs1680636 and CXCL12 in cardiac fibroblasts, overlaid with GRiNCH TAD calls. [Images not available. See PDF.]

Diamonds in Hi-C contact maps represent significant interactions (FDR < 0.1) between binned regions of the genome based on the HiCCUPS algorithm. The HiCCUPS statistical test corresponds to a one-sided Poisson test comparing observed vs. expected contacts in four local neighborhoods (donut, horizontal, vertical, and lower-left). BH-FDR multiple testing correction is applied within lambda-bins. Dotted lines in Manhattan plots indicate genome-wide statistical significance thresholds of 5E-8 for one-sided Chi-square tests of association between SNP genotypes and phenotypes. TADs topologically associating domains; ATAC assay for transposase-accessible chromatin; SNP single nucleotide polymorphism; LD linkage disequilibrium; CEU Utah residents (Center d′Etudes du Polymorphisme Humain) with Northern and Western European ancestry; GWAS genome-wide association study.

rs7549250 and rs4129267 both map to intronic regions of cytokine receptor IL6R, and their proxy SNPs have been associated with multiple heart disease-related traits. Proxies for rs7549250 were identified in GWASs of myocardial infarction (MI) (rs12118721, OR 1.06 [1.04-1.08]⁴⁰, P = 1E-07) and coronary artery disease (CAD) (rs6689306, OR 1.06 [1.04-1.08], P = 3E-09⁴⁰; rs6689306, OR 1.05 [1.03-1.07], P = 2E-09⁴¹; rs4845625, OR 1.0458 [1.0453-1.0463], P = 7E-15³⁸), while a proxy SNP of rs4129267 was associated with decreased C-reactive protein levels (rs4537545, -11.5% [-14.4% to -8.5%, P = 1.3E-12) (Table 1, Fig. 5). Both rs7549250 and rs4129267 overlap with several transcription factor footprints, and interact with IL6R TSSs ~20-kb away in cardiac fibroblasts. While there were no large GRiNCH TAD calls at this locus, other algorithms identified large-scale TADs nearby (Figure S5). rs7549250 is an IL6R eQTL in the left ventricle (P = 1.6E-20) and atrial appendage (P = 1.4E-06), while rs4129267 does not tag any IL6R heart eQTLs but is a significant eQTL of IL6R in tibial artery (P = 2.4E-09), whole blood (P = 4.6E-07), esophagus muscularis (P = 1.9E-06), and transverse colon (P = 1.5E-05). However, IL6R is lowly expressed in both primary and immortalized fibroblasts.

Fig. 5 Significant interactions between myocardial infarction- and coronary artery disease-associated open chromatin at rs7549250 and rs4129267 with IL6R in cardiac fibroblasts, overlaid with GRiNCH TAD calls. [Images not available. See PDF.]

Diamonds in Hi-C contact maps represent significant interactions (FDR < 0.1) between binned regions of the genome based on the HiCCUPS algorithm. Dotted lines in Manhattan plots indicate genome-wide statistical significance thresholds of 5E-8 for one-sided Chi-square tests of association between SNP genotypes and phenotypes. TADs topologically associating domains; ATAC assay for transposase-accessible chromatin; SNP single nucleotide polymorphism; LD linkage disequilibrium; CEU Utah residents (Center d′Etudes du Polymorphisme Humain) with Northern and Western European ancestry; GWAS genome-wide association study.

rs35346340 is in LD with two heart disease GWAS SNPs: rs2521501 was associated with CAD (OR = 1.06 [1.04-1.08, P = 5.0E-08) and MI (OR = 1.07 [1.04-1.09], P = 1.5E-07)⁴⁰, while rs4932373 was associated with MI only (OR = 0.93 [0.92-0.94])³⁸ (Table 1, and Fig. 6). rs35346340 also resides in cardiac fibroblast ATAC-seq footprints corresponding to the binding sites of multiple TFs, including KLF15. Interestingly, although rs35346340 is in the first intron of FES, a gene adjacent to FURIN, this SNP only interacts with the FURIN in our Hi-C data. GRiNCH, directionality index, and insulation score did not find any TADs at this locus, but the arrowhead algorithm found larger TADs in this region (Figure S6). rs35346340 is a significant FURIN eQTL in esophageal mucosa (P = 2.0E-13) and a nominally significant eQTL of this gene in aortic artery (P = 1.0E-04), and FURIN is highly expressed in both primary and immortalized cardiac fibroblasts.

Fig. 6 Significant interaction between myocardial infarction- and coronary artery disease-associated open chromatin at rs35346340 with FURIN in cardiac fibroblasts, overlaid with GRiNCH TAD Calls. [Images not available. See PDF.]

Diamonds in Hi-C contact maps represent significant interactions (FDR < 0.1) between binned regions of the genome based on the HiCCUPS algorithm. The HiCCUPS statistical test corresponds to a one-sided Poisson test comparing observed vs. expected contacts in four local neighborhoods (donut, horizontal, vertical, and lower-left). BH-FDR multiple testing correction is applied within lambda-bins. Dotted lines in Manhattan plots indicate genome-wide statistical significance thresholds of 5E-8 for one-sided Chi-square tests of association between SNP genotypes and phenotypes. TADs: topologically associating domains; ATAC: assay for transposase-accessible chromatin; SNP: single nucleotide polymorphism; LD: linkage disequilibrium; CEU: Utah residents (Center d′Etudes du Polymorphisme Humain) with Northern and Western European ancestry; GWAS: genome-wide association study.

Considering the importance of cardiac muscle to heart disease, we extended integrative Hi-C, ATAC-seq, and bulk RNA-seq analyses to commercially available primary human cardiac myocyte (HCM) (Tables S1-7) Since terminally differentiated HCM cells no longer proliferate and are experimentally challenging to culture, we instead characterized progenitor, non-contractile HCM cells, which can be propagated in-vitro and still express cardiomyocyte marker genes (SLC8A1, PLN, TTN) and transcription factors (MEF2A, MEF2C) more strongly than cardiac fibroblasts (Figure S1B, C). In total, 3982 active and functional SNP-gene pairs were identified in HCMs, which link 1,498 unique accessible SNPs to promoter regions (Figure S7A, Supplementary Data 2). Of these SNP regions, 137 SNPs were in binding sites of expressed transcription factors (Figure S7B). Of all the SNP-gene pairs identified in HCM and NHCFV, 1006 SNP-gene pairs were active in both cell types and 750 target genes were shared, suggesting that certain disease-associated regulatory links may be active in multiple cell types (Figure S7C). Notably, in HCM cells we also identified the rs1680636-CXCL12 and rs35346340-FURIN SNP-gene interactions we observed in NHCFV (Figure S8).

We further confirmed that identified target genes were expressed in diseased adult human cardiac fibroblasts by querying human cardiac myopathy single-nuclear RNA-seq (snRNA-seq) and spatial transcriptomic datasets, and explored how expression changes across cell types, disease states, and zones in the heart^42,43. All five target genes and corresponding heart disease-associated transcription factors involved in the prioritized SNP-gene interactions were expressed in both fibroblasts (3–33% of cells for target genes, 1–36% of cells for transcription factors) and cardiomyocytes (1–26% of cells for target genes, 0.1–35% of cells for transcription factors) (Figure S9A). Of the target genes, fibroblasts more highly express CXCL12, TBC1D32, IL6R and FURIN, while cardiomyocytes are the highest expressors of GJA1 of any cell type in the heart. Comparison of target gene expression in snRNA-seq and spatial data revealed differing gene expression levels among acute myocardial infarction (AMI) patients compared to ischemic and non-cardiomyopathy patients (Figure S9B, C). In addition, we observed co-localization between fibroblasts/cardiomyocytes and each target gene in spatial transcriptomic human myocardial infarction data (Figure S9D, E). Expression is generally concentrated in the ischemic zone of AMI patients and less prevalent in control patients (Figure S9F). Together, the variation in gene expression across disease types and heart regions suggests that tissue-dependent contexts, in addition to genetic variation, coordinate to regulate target gene expression in cardiac fibroblasts and myocytes.

Validation of SNP-function relationships using single-cell profiling of CRISPR-mediated enhancer deletions

To measure the effect of non-coding, disease-associated loci on predicted target genes as well as their broader transcriptomic consequences, we used a modified direct capture Perturb-seq approach to knock out genomic loci in immortalized human cardiac fibroblasts (iHCF) (Fig. 7A). For each locus, two guide RNAs flanking each side of the predicted TF binding footprint were designed to induce deletion of the TF footprint. Upon recruitment of Cas9 to the loci, CRISPR-mediated deletion would excise and inactivate the TF footprint region and, consequently, modulate genes that are regulated by these regions, which are located between 16 to 851 kilobases away from SNPs (Figure S10).

Fig. 7 Validation of enhancer regulation of target genes in human cardiac fibroblasts using single-cell Perturb-seq. [Images not available. See PDF.]

A Diagram of gene perturbation strategy to delete enhancer regions and perform Perturb-seq. DSB: double-stranded DNA break. B Stacked bar plot showing percentage of cells with effective KO alteration as inferred by Mixscape. NP non-perturbed; KO knockout. C Linear-discriminant analysis UMAP projections of cardiac fibroblasts, labeled by gene perturbation (left) or by cell cycle state (right). D Violin plots showing gene expression changes in predicted target genes of specific enhancer sites (showing 1000 cells per group). Exon-targeting ACTA2 gRNA was used as a positive control for gene editing activity. P values were calculated using the two-part likelihood ratio test employed by MAST differential expression analysis. E Scatter plots showing transcriptome-wide changes in per-cell detection rate and log-fold-change between a perturbation class and the non-targeting control (NTC) group. Black gene labels show the top differentially expressed genes that are unique to each perturbation. Red gene labels show predicted target genes of enhancer regions. F Gene module analyses of genetically edited fibroblasts using gene markers of fibroblast cell states during human heart failure (myofibroblast differentiation, APOE⁺ fibroblasts, inflammation, and fibrosis). 1000 cells were subsampled per group. P values were calculated using the two-sided Wilcoxon rank sum test. *P < 0.05, **P < 0.001, ***P < 0.00001, n.s. not significant. Elements in Fig. 7A were created using BioRender. Lu, Daniel R. (https://BioRender.com/rzk0wrj).

Since primary cardiac fibroblasts are unable to continuously proliferate after twenty passages, we used low-passage iHCFs for Perturb-seq, which can maintain cell doubling times after Cas9 selection and gRNA selection. Separate iHCF cultures were transduced with a pool of guides against putative enhancer regions (ER) at chromosome regions 1q21.3, 6q22.31, 10q11.21, and 15q26.1. We also included a separate pool of guides targeting exonic regions of ACTA2, which is highly expressed in cardiac fibroblasts and served as a positive control. Twenty-three days after gRNA and Cas9 co-transduction, cell cultures were harvested and processed to measure transcriptome-wide mRNA expression and guide RNA expression.

After filtering cell transcriptomes for gene quality, we detected 7519 to 16,436 cells for each perturbation, across two replicate libraries per perturbation group. We detected a high abundance of gRNA expression within cells transduced with gRNA (average of 200-900 UMI per cell) but not with cells only receiving the gene-editing enzyme. Average guide RNA expression levels per cell were typically 70–200 times higher than endogenous mRNA levels, and > 95% perturbation efficiency was observed for all groups (Fig. 7B). After regressing out cell cycle effects (Fig. 7C) and clustering cells by perturbation group, we performed single-cell differential expression analysis. Cells receiving ACTA2 gRNA showed a significant decrease in ACTA2 mRNA expression (Fig. 7D, E), confirming effective gRNA delivery and Cas9 activity. For all groups of cells receiving guides, including the non-targeting guides, we observed differential expression of DNA-sensing IFN and OAS family genes, and hence excluded these genes from consideration.

For regions ER6q22.31, ER10q11.21, and ER15q26.1 whose predicted target genes were GJA1/TBC1D32, CXCL12, and FURIN, respectively, we observed significant differential expression of at least one target gene following gene editing at each predicted regulatory element (Fig. 7E). CXCL12 and FURIN were significantly downregulated following deletion of ER10q11.21, and ER15q26.1, while GJA1 was upregulated after we deleted ER6q22.31. When comparing NTC vs Cas9-only cells, these target genes were either not significantly differentially expressed or showed opposing differential expression trends, suggesting that expression changes observed following gene editing were a direct result of regulatory element deletion, rather than artifacts of cell stress associated with cell transduction and antibiotic selection. Excision of region ER1q21.3 did not significantly downregulate predicted target gene IL6R, nor did we detect changes in TBC1D32 expression after ER10q11.21 knockout. The extremely low expression levels of both IL6R and TBC1D32 from bulk (Table 1, Figure S11A-B) and single-cell RNA-seq (Figure S7) indicate that the specific cell line used may not be a reliable model system to interrogate regulation of these particular genes.

To ensure that disruption of enhancer sites did not alter gene expression due to non-specific alterations in chromatin structure, we examined whether genes most proximal to the targeted perturbation sites that were not predicted regulatory targets were differentially expressed. For all sites, the most proximal genes (SHE and TDRD10 in ER1q21.3; HSF2, SERINC1, and PKIB in ERq22.31; TMEM72 and ZNF485 in ER10q11.21; and MAN2A2 in ER15q26.1) were not differentially expressed. Also, none of the predicted target genes were differentially expressed in our ACTA2 and NTC control perturbations. Together, these findings suggest that the transcriptional consequences we observe were caused by enhancer disruption and not non-specific changes in chromatin structure.

KEGG-based pathway analysis of differentially expressed genes following knockout found enrichment of additional genes and pathways relevant to heart failure (Table S7, and Figure S11C). Following deletion of enhancer region ER10q11.21 (predicted to interact with GJA1 and TBC1D32), pathway analysis identified COX7A1 downregulation in multiple heart disease-relevant pathways (cardiac muscle contraction (hsa04260), diabetic cardiomyopathy (hsa05415), and vascular smooth muscle contraction (hsa04270)). In addition, deletion of the enhancer interacting with FURIN (ER15q26.1) resulted in decreased expression of its substrate, MMP2, and of ACTG1, a core enrichment gene in dilated (hsa05414) and hypertrophic (hsa05410) cardiomyopathy pathways. No heart disease-related KEGG pathways were detected from the differentially expressed genes of ER1q21.3 and ER10q11.21 knockout cells.

We next performed in-depth analyses of transcriptional changes induced by CRISPR KO of predicted enhancers at loci associated with coronary disease. Using multi-gene signatures previously defined by single-cell RNA-seq of fibroblast cell states in acute myocardial infarction^42,44, we observed changes in genes associated with fibrotic (POSTN, COMP, COL1A1, THBS4, COL3A1, FBN1, PPRX, FOSX1, MEOX1, RUNX1, EDNRA), inflammatory (CCL2, CCL11, THBS1, PTGDS, GPC3), and myofibroblast (ACTA2, TAGLN) states after CRISPR KO (Fig. 7F). Of note, knockout of ER15q26.1, which reduced expression of myofibroblast-promoting protease FURIN, also reduced expression of myofibroblast gene signatures. We also observed changes in APOE⁺ fibroblast signatures due to perturbation, which have been phenotypically linked to cholesterol dysregulation and fibrosis⁴⁵. We asked whether changes in cholesterol regulation can be inferred by applying cholesterol-related gene signatures to perturbed cell states and found perturbations that enrich for the APOE⁺ fibroblast signature (i.e. in ER1q21.3) associated with stronger differential regulation of cholesterol metabolic and homeostatic processes compared to perturbations that reduced the APOE⁺ fibroblast signature, suggesting a possible role for ER1q21.3 regulation on APOE expression and cholesterol metabolism (Figure S9D). These results suggest that perturbation of these sites has the potential to affect the expression of interacting genes (identified by Hi-C) as well as downstream effectors that may alter cardiac fibroblast function.

Overall, our Perturb-seq screen validates the regulatory relationships between three of four predicted enhancer regions harboring heart disease-associated GWAS loci: 6q22.31-GJA1 (atrial fibrillation), 10q11.21-CXCL12 (coronary artery disease), and 15q26.1-FURIN (myocardial infarction and coronary artery disease). We also identified possible roles for these enhancer sites in modulating global gene expression programs associated with fibrosis, inflammation, and cholesterol regulation, reinforcing the direct pathophysiological roles of these non-coding sites on cellular phenotypes.

Discussion

In this study, we traced regulatory circuits from heart disease SNPs to their target genes, by integrating genome-wide interaction, epigenomic, and transcriptomic data from human cardiac fibroblasts. We identified 3D chromatin interactions using both GOTHiC and HiCCUPs to ensure that we captured both short- and long-range loops⁴⁶. Notably, across both Hi-C loop calling tools, over 2/3 of SNP-gene interactions prioritized for regulatory relationships involve target genes that were not mapped to the original SNP, a proportion consistent with findings from orthogonal approaches⁴⁷. ATAC-seq provided evidence to suggest that GWAS loci reside within putative regulatory elements, and transcriptomic profiling confirmed the cell-specific expression patterns of target genes. Deletion of predicted enhancer regions demonstrated that these regions actively regulate cellular transcriptional programs and, in some cases, are sufficient to modulate target gene expression. Our multi-omic, integrative approach revealed disease genes that otherwise would not have been identified by considering the genes most proximal to GWAS signals alone.

Our target genes and the TFs corresponding to binding motifs in ATAC-seq footprints were enriched for several known cardiac remodeling and fibrosis genes and TFs, further affirming the value of integrating GWAS signals with multi-omics data from relevant tissues. Such integrative multi-omics approaches are especially important for heart failure, where disease heterogeneity in addition to other factors have precluded a comprehensive understanding of its genetic basis, particularly polygenic forms of heart failure⁴⁸. Below, we summarize the evidence suggesting pathophysiological relevance of several of our identified transcription factors and target genes.

Transcription factor FOXJ3’s Target Genes GJA1 and TBC1D32

Gap Junction Protein Alpha 1 (GJA1) encodes Connexin 43 (Cx43), a gap junction component involved in pathogenic trans-differentiation of fibroblasts into myofibroblasts, which migrate to injured areas to deposit extracellular matrix. Specifically, Cx43 mediates Transforming Growth Factor B (TFG-β)-induced fibroblast-to-myofibroblast phenoconversion in rats, and Cx43 knockdown and overexpression respectively inhibited and stimulated the promoter activity of α-SMA (a myofibroblast marker) in cardiac fibroblasts⁴⁹. Furthermore, another study found that knockdown of Connective Tissue Growth Factor (CCN2/CTGF), which encodes a primary target and downstream effector of TGF-β, reduced expression of Cx43⁵⁰. These prior findings and our own prioritization results are all consistent with the notion that GJA1 plays a key role in heart failure-associated myocardial fibrosis. TBC1 Domain Family Member 32 (Tbc1d32) was identified in a recessive forward genetic screen in fetal mice in which the authors performed dense phenotyping for congenital heart defects and exome sequencing to identify recessive mutations. Specifically, the Tbc1d32 mice showed laterality defects that may be related to sonic hedgehog (SHH) signaling transduced or modulated by the cilium⁵¹. Although deletion of ER_6q22.31 did not alter TBC1D32 expression, we detected decreased expression of COX7A1, which encodes the contractile muscle-specific and major cardiac isoform of cytochrome c oxidase (Cox). Homozygous and heterozygous Cox7a1 knockout mice demonstrate dilated cardiomyopathy⁵². While these target genes will require deeper characterization to clarify their therapeutic potential, our analyses point to the key gene targets, TFs, and pathways implicated at this atrial fibrillation-associated GWAS locus.

Target gene CXCL12 and transcription factor ATF3

C-X-C motif chemokine ligand 12 (CXCL12), also known as Stromal Cell Derived Factor 1 (SDF1), is a ligand for C-X-C motif receptor 4 (CXCR4), a G protein-coupled receptor with roles in angiogenesis, myocardial ischemia, and injury-induced restenosis⁵³. In our data CXCL12 interacts with CAD-associated rs1680636, which is in a binding motif of Activating Transcription Factor 3 (ATF3), a basic leucine zipper family TF involved in cardiac remodeling. Specifically, ectopic expression of Atf3 in Atf3 knockout mice inhibited angiotensin II-induced cardiac fibrosis and hypertrophy⁵⁴. Furthermore, previous cancer studies have implicated an interaction between Atf3 and CXCL12 in tumor progression^55,56. Taken together, these results suggest CXCL12 may be an important target of Atf3 involved in the etiology of heart failure and associated cardiovascular outcomes.

Target gene IL6R and transcription factor KLF12

The Interleukin 6 Receptor (IL6R) ligand IL-6 is a cytokine involved in cell growth and differentiation, as well as the immune response, and is tightly associated with the liver-derived inflammatory biomarkers C-reactive protein (CRP) and fibrinogen. All three are significantly associated with increased risk of heart disease⁵⁷. IL-6 may mediate cardiac remodeling consequent to Ang II, as IL-6 knockout mice were resistant to Ang II-induced cardiac dysfunction, myocardial inflammation, and fibrosis compared to wild-type mice⁵⁸. Consistent with the notion that IL6R may be directly involved in inflammation that leads to heart disease, a meta-analysis found that a non-synonymous SNP in IL6R (rs8192284 [p.Asp358Ala]) and its proxy rs7529229 are robustly associated with coronary heart disease (per-allele OR 0.95 [0.93-0.97], P = 1.5E-05), as well as increased IL-6 and decreased levels of downstream effectors CRP and fibrinogen⁵⁹. Interestingly, another significant SNP at this locus, rs7549250, is in a TF footprint of Kruppel-like factor 12 (KLF12), which is associated with rheumatoid arthritis. Tocilizumab, which is commonly prescribed for rheumatoid arthritis and other chronic inflammatory diseases, competitively inhibits IL6R, and, like rs8192284, is associated with increased IL-6 and decreased CRP and fibrinogen. However, the drug has a directionally inconsistent effect on cardiovascular disease due to its pro-atherogenic lipid profile⁵⁹. These findings suggest crucial roles for IL6R in heart disease and point to a need to more deeply characterize KLF12 targets and mechanism(s) of action in cardiac-relevant tissues.

Transcription factors EGR1, KLF15, and SNAI1, and target gene FURIN

Another SNP, rs35346340, is in the TF footprints of several TFs which have been associated with phenotypes related to heart failure: Early Growth Response 1 (EGR1), Kruppel-like factor 15 (KLF15), and Snail 1 (SNAI1). EGR1 is a C2H2-type Zn finger transcription factor and Egr1 knockout mice show increased cardiac fibrosis⁶⁰. In cardiac fibroblasts, KLF15 inhibits expression of CCN2/CTGF⁷, thereby preventing cardiac fibrosis. SNAI1 is a Zn finger transcriptional repressor, and Snai1 was co-expressed with the fibrosis marker periostin in post-MI mice, suggesting a role in post-MI de novo fibrosis⁶¹. The target gene of rs35346340, Furin, Paired Basic Amino Acid Cleaving Enzyme (FURIN), is important to the renin-angiotensin system and sodium-electrolyte balance, and was prioritized in an integrative polygenic analysis of GTEx aorta eQTLs and GWAS data as the likely causal gene underlying a systolic blood pressure association signal⁶². Expression of FURIN has also been found to be elevated in canines⁶³ and mice⁶⁴ with heart failure. The lower FURIN expression and reduced inflammatory signature we observed after rs35346340 knockout suggest that the enhancer at this locus promotes a fibroblast phenotype that is detrimental during heart failure. Intriguingly, we found that deletion of the enhancer interacting with FURIN also resulted in differential expression of γ-actin (ACTG1), which is the major actin isoform in the heart and important in cardiac muscle contraction⁶⁵.

As mentioned, our Hi-C data, even at high 2-kb resolution, could not always resolve to a single gene’s TSS, depending on local genetic architecture. Hi-C is also limited to assaying pairwise interactions, and we thus may miss further characterizations of genes and variants with multiple interaction effects (although some proportion of these will be found in orthogonal pairwise interactions) (Supplemental Notes). Additionally, our use of ATAC-seq data to identify enhancer regions is by no means comprehensive. Some enhancers are only active in stimulated cells, as evidenced by studies of immune cells^66,67, and may not have been represented in our ATAC-seq data of wild-type primary cells. Furthermore, not all intergenic ATAC-seq-derived open chromatin peaks necessarily represent enhancers, with some likely harboring other types of cis regulatory elements with different impacts on gene expression. A limitation of our validation work is the use of immortalized cells over a primary line which are deemed to be more biologically relevant. Numerous prior studies have utilized immortalized cell lines for CRISPR based screening given the practical advantages of scalability^68,69. Lastly, addition of other types of epigenomic (i.e. histone modification) data could aid in more deeply characterizing the direction and nature of regulatory relationships underlying GWAS signals.

Methods

Ethics

The design and conduct of this study complies with all relevant regulations regarding the use of human tissues and data, in accordance to the criteria set by the Declaration of Helsinki. All genetic association data analyzed were anonymized population-level summary statistics and therefore did not contain any personal identifying information. All cell lines were commercially obtained from normal non-diseased adults. Ethnicity and age information was not available from vendors due to privacy concerns or incomplete records. For the above reasons we did not require institutional review/ethics board approval for this study.

Cells

Primary human fibroblasts and cardiac myocytes

We generated Hi-C, ATAC-seq, and RNA-seq data from primary cells to model heart failure-relevant cellular contexts. Specifically, we cultured 50 million primary Normal Human Cardiac Fibroblasts—Ventricular (NHCFV) as well as Human Cardiac Myocytes (HCM) primary cells, which are commercially available from Lonza (Catalog Number CC-2904) and PromoCell (Catalog Number C-12810) respectively.

Immortalized human cardiac fibroblasts

We performed gene editing experiments using Immortalized Human Cardiac Fibroblasts generated from adult heart ventricle primary cells (iHCF), which are commercially available from Applied Biological Materials Inc (Catalog Number T0446). These experiments were performed in iHCF as opposed to primary NHCFV because compared to cell lines, primary cells have limited life span and are less tolerant of manipulations such as electroporation and gene editing. To account for this difference and enable comparisons of peak calls and gene expression levels between immortalized and primary fibroblasts, we also generated ATAC-seq and RNA-seq data from this cell line.

NCHFV, HCM, and iHCF cells were obtained from male donors, and we confirmed their sex based on the presence of Y-chromosome peaks in ATAC-seq data from each sample. The cell lines’ respective manufacturers tested for cell morphology, adherence rate, cell viability, and cell-type specific markers using flow cytometry. All cell lines tested negative for mycoplasma by manufacturers. We further tested for morphology and myocyte contractility as well as cardiomyocyte expression profiles to authenticate these cells.

Hi-C library construction and data processing

We contracted Arima Genomics to construct in-situ Hi-C libraries for 10 million NHCFV and HCM cells based on the protocol from Rao et al. ²⁷. For NHCFV cells, they used an enzyme cutting at ^GATC, while a mixture of restriction enzymes that cut at ^GATC and G^ANTC was used for HCM cells. We then sequenced over 2 billion paired-end reads (2.4E9 for NHCFV and 2.7E9 for HCM) using an Illumina HiSeq 4000. Reads from both biological replicates were processed together for each cell type (aligned to hg38 and read- and fragment-filtered to identify valid interaction pairs) using both a) the HiC-Pro pipeline version 2.8.1⁷⁰ and b) the Juicer pipeline version 1.5.6⁷¹ on the Juicer Amazon Machine Image (AMI) on Amazon Web Services (AWS). We called genome-wide loops representing predicted genome-wide pairwise physical interactions between 2-kb regions of the genome using GOTHiC version 1.12.0⁷² and HiCCUPS included with Juicer 1.5.6⁷¹.

Hi-C loop & TAD calling

GOTHiC

The BAM file from HiC-Pro containing valid interaction pairs (filtered for invalid interaction products), was converted to GOTHiC format using hicup2gothic, and split by chromosome. We ran GOTHiC in base R 3.4.4 on AWS to call significant interactions from each chromosome-specific file. Statistical significance was set at FDR < 0.05 and contact counts > 10 based on previous studies^46,73.

HiCCUPS

Using the.hic file generated by Juicer, we ran HiCCUPS on an AWS GPU instance to identify loops at multiple resolutions (2-kb, 5-kb, 10-kb, and 25-kb) by calling peaks enriched compared to four local backgrounds: (i) the donut shape around a pixel, (ii) in the lower-left quadrant of the donut neighborhood, as well as the (iii) vertical and (iv) horizontal neighborhoods²⁷. At each resolution we used the following arguments for the peak width (p), window (i), FDR threshold (f), and centroid distance (d) parameters suggested by the authors: 2-kb (p = 10, i = 20, f = 0.1, and d = 20,000), 5-kb (p = 4, i = 7, f = 0.1, and d = 20,000), 10-kb (p = 2, i = 5, f = 0.1, and d = 20,000), and 25-kb (p = 1, i = 3, f = 0.1, and d = 50,000). Vanilla coverage was used for normalization, as the Knight-Ruiz method did not converge for some chromosomes at 2-kb and 5-kb.

TAD calling

We used the .hic files generated by Juicer to infer TADs with the arrowhead algorithm at resolutions of 5, 10, 25, 50, 100, and 500 kb, using Vanilla coverage normalization. To ensure robust results, we also converted our .hic files to sparse format at these same resolutions, and subsequently applied GRiNCH⁷⁴, with a neighborhood size of 5 bins and default parameters, to infer TADs via graph regularized non-negative matrix factorization and clustering. TADtool⁷⁵ was also used to infer TADs using both directionality index and insulation score. We discovered TADs using a window size 5x the resolution, and a cutoff value per-chromosome representing the top 5% of the observed or absolute values for insulation score and directionality index, respectively. The TAD intervals identified in each cell type were then overlapped with the start and end of our identified SNP-gene pairs using bedtools⁷⁶ version 2.22.1. We utilized bedtools intersect with parameter -f 1.0 to assess how many SNP-gene pair bin intervals fell entirely within the same TAD.

Bulk ATAC-seq data processing, peak calling, transcription factor binding site enrichment and footprinting

We constructed single libraries from 1 million NHCFV cells and 100,000 iHCF cells, and two ATAC-seq libraries from 1 million HCM cells, using the protocol from Buenrostro et al.⁷⁷. Libraries were sequenced at > 150 million paired-end reads per sample on either an Illumina HiSeq 4000 or NextSeq 500, resulting in > 50 million de-duplicated retained reads per sample. We used the minimum number of iHCF cells required for ATAC-seq (10x fewer than the primary cells) due to the difficulty of culturing these cells, and to save material for other experiments. We processed raw reads and called open chromatin regions using the ENCODE ATAC-seq pipeline⁷⁸, setting hg38 as the reference genome. We retained putative open chromatin regions from naïve overlap peaks called from two pseudo-replicates for downstream analyses.

To find cell-specific transcription factor binding sites (TFBSs) in ATAC-seq peaks, we performed a nucleotide-level computational search for transcription factor (TF) footprints using HMM-based IdeNtification of Trancription factor footprints (HINT)-ATAC from Regulatory Genomics Toolkit (RGT) release 0.11.4^{79, 80, 81–82}. The output included motif-predicted binding sites of known TFs from the HOCOMOCO⁸³ and JASPAR^84,85 databases, as well as bit-scores of motif matching^80,86. TF motif enrichment analysis was performed using Hypergeometric Optimization of Motif Enrichment (HOMER)⁸⁷. We also analyzed differentially footprinted regions to compare predicted TF binding between cell types.

RNA sequencing and transcript quantification

Global transcript expression in NHCFV, iHCF, and HCM cells was assessed by RNA sequencing (RNA-seq). RNA-seq was performed on a cDNA library prepared from total RNA (2 μg; RIN score > 9.5) of 100,000 cells isolated using Mirvana miRNA RNA isolation kits (Ambion, Grand Island, NY) with on-column DNase treatment. Total RNA quality and concentration were determined using the Bioanalyzer (Agilent, Santa Clara, CA) and Nanodrop (ThermoScientific, Wilmington, DE). cDNA was prepared using a modified protocol based on the Illumina Truseq RNA Sample Preparation Kit (Illumina, San Diego, CA) and the published methods for strand-specific RNA-Seq^88,89. After size selection of libraries (Pippen Prep; SAGE Biosciences, Beverly, MA), dUTP-containing cDNA strands were destroyed by digestion of USER enzymes (New England Biolabs, Ipswich, MA) followed by PCR enrichment for introduction of strand specificity. The enriched cDNA libraries were analyzed in Agilent Bioanalyzer and quantified by Quant-iTTM Pico-Green assays (Life Technologies). RNA sequencing reads (Illumina HiSeq platform, 75 bp paired-end sequencing) were aligned to hg38 and Fragments per Kilobase per Million sequenced, Quantile normalized (FPKQ) values were determined using Array Suite software (Omicsoft, Cary, NC) and in-house software. Genes with FPKQ values > 0 were considered expressed.

Genome annotations

NHGRI SNPs

We queried the NHGRI-EBI GWAS Catalog⁹⁰ using the search terms “heart disease,” “obesity,” and “pulmonary hypertension” on September 5, 2019, and also included a heart failure GWAS preprint that is now published⁹¹. For each SNP, we noted the proximal gene(s) based on the annotations “Mapped Gene” (the gene(s) mapping to intragenic SNPs, or genes up- and down-stream of intergenic SNPs) and “Reported Gene” (gene(s) reported by the study author, typically based on proximity to the reported SNP). We then used LDproxy from LDlink release 3.7.2⁹² to identify all proxy SNPs (r² > 0.8 in the European populations).

FANTOM5 TSSs

To exhaustively identify all human transcription start sites (TSSs), we downloaded Biomart annotated hg38 CAGE Peaks from across all human cell types and tissues assayed by the FANTOM5 consortium⁴⁵, and intersected these TSSs with gene symbols from HGNC⁹³. In total we obtained 96,254 unique TSSs mapping to 20,224 unique genes.

GTEx eQTL data

We downloaded single-tissue cis-eQTLs computed from 44 human tissues in GTEx release v8⁹⁴. We first annotated eQTL SNPs with LD-based loci to cluster them into locus-gene pairs, which we then used to annotate our Hi-C-derived SNP-gene pairs in the same loci as having eQTL support.

Data integration and analysis

Our analytical approach aimed to identify and use potentially functional SNPs underlying disease GWAS signals, layered with functional genomic information (Hi-C, ATAC-seq, and RNA-seq), to trace regulatory circuits between putative enhancers and their target genes.

1. Due to linkage disequilibrium (LD), GWAS loci typically contain multiple significant disease-associated SNPs indistinguishable by effect-size and p value, and it is unclear which of these SNPs (or proxy SNPs not genotyped in the original GWAS) underlie the association signals. We therefore studied all reported GWAS SNPs, as well as their proxies (r² > 0.8 in the CEU population). To prioritize potentially functional SNPs based on overlap with putative enhancers, we layered accessible chromatin regions and noted footprints containing motif-predicted transcription factor (TF) binding sites, especially those of expressed TFs.

2. We used 2-kb resolution Hi-C data to identify significant pairwise interactions between the SNPs in predicted enhancers from (i) and accessible transcription start sites (TSSs). We hypothesized that SNPs disrupting the function of TSS-interacting enhancers may underlie their significant disease associations.

3. SNPs in the SNP-TSS pairs from (ii) that are (or are in LD with) significant GTEx eQTLs of the target gene have additional evidence to suggest that the SNP may affect expression of the TSS’ target gene. However, GTEx eQTLs were generated from whole-tissue mixtures of cells and may not always be concordant with SNP-TSS pairs identified from the above analyses of a homogeneous primary cell population.

This approach identifies putative enhancers whose perturbation by disease-associated common genetic variations may affect regulation of genes involved in the etiology of complex diseases. Further, this method distinguishes potential causal SNPs from their proxies at disease-associated loci.

Data visualization

We visualized physical interactions between accessible SNPs and TSSs using Juicebox version 2.15⁹⁵ and Hi-Glass version 2.5.7⁹⁶. We plotted Hi-C and GWAS data using PlotGardener⁹⁷. Hi-C heatmaps depict observed/expected interaction counts, to accentuate distal interactions and reduce noise from proximal interactions along the x-axis. Since small-scale TADs ( < 50-kb) appeared redundant with loops, we only plotted larger TADs (50-, 100-, and 500-kb) to improve the clarity of our figures (Supplemental Notes).

Experimental validation of high-priority SNP-gene pairs

Generation of dual-flanking guide RNA libraries for Perturb-seq

We used the Perturb-seq gRNA-mediated KO strategy in the ventricular cardiac fibroblast cell line (iHCF) to introduce deletions at a subset of ATAC-seq regions prioritized for their overlap with heart disease-related TF footprints and interactions with accessible TSSs of heart disease-related genes. To improve deletion efficiency, we introduced two or more single-stranded guide RNAs (gRNAs) with homology to sequences proximal to the regions 5’ (upstream) and 3’ (downstream) of putative enhancers ranging from 746 to 2863 base pairs. Valid gRNAs were predicted using the Sigma Design Tool⁹⁸ and ChopChop algorithm⁹⁹. The gRNAs with the (1) highest predicted efficiency and specificity, (2) lowest predicted self-complementarity and promiscuity (no sites in genome that are < 3 nucleotide mismatch to gRNA), and (3) target sequences most proximal to the transcription factor binding site (0 to 580 bp) were selected for cloning into the LV15 vector (Millipore Sigma). This vector expresses gRNAs under control of the human U6 promoter and gRNA scaffolds contain Capture Sequence 2 in the stem loop, which enables hybridization onto 10X Genomics 3’ gel beads and detection via Perturb-seq. OffSpotter¹⁰⁰ and cas-offinder^101,102 were used to check gRNA sequences for off-target activity.

The gRNAs were designed to flank each predicted transcription factor binding site such that guides flanking the 5’-end were inserted into an LV15 vector backbone (Millipore Sigma) containing a Puromycin-resistance-BFP (blue fluorescent protein) selection marker, and guides flanking the 3’-end were inserted into an LV15 vector backbone containing a Zeocin-resistance-BFP selection marker. In total, four gRNAs were utilized for each enhancer locus. gRNAs against the exon regions of ACTA2 were also designed as a positive control for targeted deletion, and two non-targeting gRNAs were selected from a published library¹⁰³ as a negative control.

Lentivirus production and quality control

Cloned plasmids containing gRNA protospacer sequences were packaged into lentivirus at a minimum formulation of 10⁷ viral particles per mL in 200 uL (Millipore Sigma). Plasmid quality and yield were initially verified by plasmid restriction digest and p24 antigen ELISA titer. A functional titer was also calculated to determine lentivirus amounts for iHCF transduction.

Guide RNA library titration for Perturb-seq

1 mL of 2 × 10⁵ cells/mL iHCF (ABM Cat# T0446) cell suspension prepared with FibroGRO^TM-LS complete media (Millipore Sigma Cat # SCMF002) + 8 µg/mL polybrene (Millipore Sigma Cat # TR-1003-G) was placed into each well of a 12-well plate. Lentivirus gRNA library was thawed at room temperature and mixed by gently tapping the tube several times with fingers. 100 µl of 10-fold dilution of the lentivirus gRNA library was prepared with FibroGRO^TM-LS complete media+ 8 µg/mL polybrene. 50 µl of 10-fold diluted lentivirus was then used to prepare 50 µl of additional 2-fold serial dilutions ranging from 1/2 to 1/1024. 25 µl of media was added to one well as negative control. 25 µl of 10-fold or additional 2-fold diluted lentivirus was added to the rest of the wells and mixed thoroughly by pipetting up and down a few times. The plate was then incubated at 37 °C for 18–20 h. The media containing virus was then removed and replaced with 1 mL of FibroGRO^TM-LS complete media. The plate was then cultured for an additional 2 days. On day 3 post lentivirus transduction, the media from each well was removed. 200 µl of trypsin was added to each well to detach the cells for 3 min at 37 °C. 800 µl of complete media was added to each well to neutralize trypsin. The cells were then transferred to 15 mL conical tubes for centrifugation at 300xg for 5 min. Each cell pellet was resuspended in 500 µl BD Cytofix^TM Fixation Buffer (DB Biosciences Cat # 554655) for 10 min at room temperature. The BFP expression of cells in each well was analyzed by FACS. We used the wells that had between 5% and 20% of cells expressing BFP to determine the actual titer. The actual titer of the lentivirus gRNA library is calculated by the following formula: titer = [(F x C_n)/V] x DF.

F: The frequency of GFP-positive cells determined by flow cytometry (0 ~ 1)

C_n: The total number of cells transduced

V: The volume of the inoculum (mL)

DF: The virus dilution factor (total inoculum volume/the volume of original lentivirus gRNA library)

After obtaining the actual titer of the gRNA lentivirus library in iHCF cells, the iHCF cells were transduced with the appropriate amount of lentivirus to achieve the desired MOI the same way as described above. The media containing virus was removed 1 day post transduction and replaced with complete media. After culturing for an additional 2 days, the cells were selected with 0.25 µg/ml puromycin for 6 days to enrich gRNA transduced cell population before harvesting cells for Perturb-seq. More than 90% of the cells were BFP positive when harvested for Perturb-seq.

Generation of CRISPR-Cas9 edited immortalized human cardiac fibroblasts

Pools of immortalized human cardiac fibroblasts (iHCF) (Applied Biological Materials, Inc.) were co-transduced with lentiviruses containing guides against the same enhancer locus and Cas9-Blasticidin lentivirus (Transomic, Inc.). Low-passage (< 5) iHCF cells were used since the cells become senescent at higher passage numbers and less able to integrate lentivirus complementary DNA. Triple antibiotic selection with Blasticidin, Puromycin, and Zeocin was applied to infected cells to ensure that cells for Perturb-seq contained the Cas9 expression cassette, ≥1 gRNA targeting the 5’-region flanking the transcription factor binding site, and ≥1 gRNA targeting the 3’-region flanking the transcription factor binding site. Guides were transduced in a semi-arrayed manner at high multiplicity-of-infection (MOI = 10) (i.e. each cell pool was transduced with two 5’ guides and two 3’ guides specific for the same locus, but the pair of 5’ and 3’ guides that infected a cell was stochastic). Our approach enabled cardiac fibroblasts to be infected at a high MOI without concern for editing of enhancer loci at multiple chromosomes in the same cell, and enabled a high fraction of cells to receive gRNA and Cas9. After applying antibiotic selection, cell viability and numbers were given time to recover over a period of 23 days, before cells were harvested for Perturb-seq. At the time of harvest, each cell pool contained greater than 100,000 live cells. Cells were harvested by trypsinization, 1x wash in culturing media, and 1x wash and resuspension into PBS + 0.04% BSA at a concentration of 300-500 cells/uL for droplet encapsulation. Viability for all cell pools was ≥ 85%.

Perturb-seq library generation and sequencing

Single-cell Perturb-seq (transcriptome and gRNA sequencing at single-cell resolution) was performed using the Chromium 3’ V3.1 method with Feature Barcoding technology for CRISPR screening (10X Genomics, User Guide CG000205 Rev D). Cells in freshly prepared PBS + 0.04% BSA were added to reverse transcription master mix with the appropriate volume of RNAse-free water per manufacturer’s guidelines. The solution containing cells and master mix was encapsulated with gel beads using partitioning oil to generate Gel Bead-in-EMulsions (GEMs) targeting 10,000 recovered cells per sample lane. After encapsulation and reverse transcription, barcoded mRNAs and gRNAs were extracted from GEMs and purified using Dynabeads (MyOne SILANE).

Next, PCR was performed using a primer cocktail (Feature cDNA Primers 1, PN #2000096) for cDNA amplification of cell-endogenous mRNA and targeted amplification of gRNAs containing CaptureSequence2. A total of 11 cycles of cDNA amplification were performed using standard temperature settings, according to manufacturer’s recommendations. During PCR clean-up using SPRIselect beads (Beckman Colter), the amplified gRNA-containing supernatant fraction was separated from the pelleted fraction, which contained amplified cDNA from endogenous mRNA. Standard library construction of 3’ gene expression libraries (GEX) of the SPRI-purified pelleted fraction was performed according to manufacturer’s recommendations, which consisted of fragmentation, A-tailing, ligation and 12 cycles of sample index PCR to append single-end 8-bp index sequences to each library pool. The gRNA-containing supernatant was SPRI-purified and amplified using two sample index PCR rounds (15 cycles and 5 cycles, respectively) to generate single-indexed final libraries for next-generation sequencing.

Completed libraries were quantified using TapeStation 4200 HSD1000 screentapes and pooled at a molar ratio of 5:1 for GEX:gRNA libraries. Libraries were sequenced on NovaSeq6000 S4 flow cells (Illumina) using the sequencing cycle specifications of 28:8:0:91 (Read1:i7 index:i5 index:Read2). We targeted a total of 75,000 combined GEX and gRNA reads per cell.

Pre-processing and filtering of Perturb-seq data

Raw sequencing BCL files were demultiplexed by GEX and gRNA library single-index (i7) barcodes to generate FASTQ files for each library. GEX FASTQ reads were aligned against human reference genome GRCh38-3.0.0, and gRNA FASTQ reads were mapped to a custom feature reference table by calling the –feature-ref argument in cellranger count (CellRanger software v6.0.0, 10X Genomics¹⁰⁴). The unfiltered UMI count matrices generated by CellRanger were used for pre-processing of data.

To assign expressed gRNA species to their respective cells, count matrices from different library samples were merged. Since each transduced cell pool only contained a subset of guides against the same target region for each target gene, we first removed all cell barcodes with non-zero expression of gRNAs that were not present in their respective cell pool. This removed 0.03%-3.8% of cell barcodes from each cell pool. We also conservatively removed cell barcodes that expressed less than 100 gRNA UMIs to remove cells with putatively low guide transcription from the analysis (the control cell pool that only expressed gene-editing enzyme was excluded from this criterion). The expression values of gRNAs were then normalized across cells using the centered log-ratio transformation and scaled¹⁰⁵. Cells were assigned to perturbation classes based on their expression of gRNAs (i.e. one class per targeted enhancer locus/gene). To remove low-quality or multiplet single-cell transcriptomes, cell barcodes that contained less than 2000 features, greater than 6500 features, or greater than 10% fraction of mitochondrial reads were discarded.

Perturb-seq analysis of transcriptional changes following CRISPR KO

Analysis of single-cell RNA-seq was performed using the Seurat package¹⁰⁶ version 3.9.9. Raw UMI count values were first transformed using log-normalization. Cell cycle phase was quantified by applying the CellCycleScoring() function to cell cycle genes identified from Tirosh et al.¹⁰⁷. Dimensionality reduction using principal component analysis (PCA) and clustering were performed to examine the confounding impact of cell cycle on transcriptional signatures following perturbation. After selecting the top 5000 most variable features using FindVariableFeatures(), expression values were scaled and the effect of cell cycle was regressed out. The first 40 principal components were used for graph-based assignment of cell neighborhoods (shared Nearest Neighbor), and clustering was performed using the Louvain algorithm (resolution=0.6) to ensure that cell clusters segregated predominantly by biological variables (e.g. gRNA class) and not by technical variables (e.g. cell cycle or library prep batch).

During gene perturbation screening, a fraction of cells may express gRNA but not undergo genetic modification. To increase the signal-to-noise ratio in Perturb-seq experiments, it is critical to distinguish cells from the same gRNA class that are functionally altered by gene editing, from cells that are not impacted by gRNA and gene-editing enzyme transduction. We used Mixscape^106,108, which uses a Gaussian mixture model to assign posterior probabilities to cells based on gRNA-mediated editing efficiency, to separate perturbed from non-perturbed cells. By specifying cells from the non-targeted control (NTC) gRNA class as a reference for a non-perturbed transcriptome, we found that ~95% of cells from our assay had a perturbed phenotype, likely due to the high MOI used in our assay. Cells predicted to be unperturbed were excluded from downstream analyses.

Differential expression for all expressed features was performed using MAST¹⁰⁹ version 1.12.0 to identify transcriptional changes induced by genomic perturbation. Due to the stringent cell filtering criteria used, we considered features with |log₂-fold-change | > 0.1 and adjusted p value < 0.05 to be statistically significant. Signature genes associated with SNP-containing enhancer regions were identified by removing any features that were differentially expressed between cells only expressing the gene-editing enzyme (editing-enzyme-expressing control) and NTC-transduced cells. We used clusterProfiler^110,111 version 4.12.0 to run gene set enrichment analysis (GSEA) to search for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that may be enriched in differentially expressed genes. Our input list of differentially expressed genes was sorted by |log₂-fold-change| in expression averaged across cells, weighted by the percentage of cells with detectable expression of that gene.

To visualize single-cell transcriptomic phenotypes, we used linear discriminant analysis (LDA) for dimension reduction to maximize variance across different gRNA classes. The gRNA class was used to distinguish cellular populations, and the first 20 components were used to calculate a UMAP projection for the cells.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Acknowledgements

We would like to thank the High Performance Computing Core Facility at the University of California, Davis, for providing computational resources during paper revision, and Ryan Potts, Simon Jackson, Scott Martin, and Bill Richards for providing helpful scientific input. We would also like to acknowledge Andy Rampersaud and Vanessa Arias for providing technical assistance. Kory Lavine would like to give thanks to the support of National Institutes of Health (R35 HL161185) and Leducq Foundation Network (#20CVD02).

Author contributions

Y.H. and R.G. designed the study and conceived the experiments. D.R.L., Y.A., J.C., C.W., J.L., T.Y., Z.J.Y., H.Z. and C.L. performed experiments. R.G., T.Y., J.M.A. and D.R.L. processed data. R.G., D.R.L., I.E. and J.M.A. analyzed data. C.L. and YH oversaw data generation and analysis. C.L., S.W., K.L., B.A. and Y.H. advised and provided resources to support the project. R.G., D.R.L. and I.E. wrote the manuscript. R.G., D.R.L., I.E., J.M.A., K.L., B.A., C.L. and Y.H. finalized the manuscript.

Peer review

Peer review information

Nature Communications thanks Richard T. Lee and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. [A peer review file is available].

Data availability

FANTOM5 CAGE Peaks are available from the FANTOM5 consortium: (https://fantom.gsc.riken.jp/5/datafiles/reprocessed/hg38_latest/extra/CAGE_peaks/). GTEx v8 bulk tissue eQTL data are available from the GTEx Portal: (https://www.gtexportal.org/home/downloads/adult-gtex/bulk_tissue_expression). ATAC-seq, Hi-C, bulk RNA-seq, and Perturb-seq data were deposited into the Gene Expression Omnibus database under accession numbers GSE281462, GSE281463, GSE281464, and GSE281465. Previously published data from Kuppe et al. (Zenodo (https://zenodo.org/records/6578047) and Amrute et al. (GSE218392) were used for analysis. We also provide source data for heart disease-related GWAS summary statistics analyzed in this study, Fig. 7, and S1, S11A–D, guide RNA oligonucleotide sequences, as well as Figure R1 showing Y chromosome ATAC-seq peaks for all primary and immortalized cells utilized in this study.

Competing interests

The authors declare the following competing interests: R.G., D.R.L., I.E., J.L., J.C., C.W., Z.J.Y., T.Y., H.Z., B.P., Y.A., S.W., B.A., C.L. and Y.H. were Amgen employees while performing experiments or analysis and own Amgen stocks. K.L. and J.M.A. were not Amgen employees and have no competing interests.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-64070-1.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Go, AS et al. Executive summary: heart disease and stroke statistics-2013 update: a report from the American Heart Association. Circulation; 2013; 127, pp. 143-152. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23283859][DOI: https://dx.doi.org/10.1161/CIR.0b013e318282ab8f]

2. Ziaeian, B; Fonarow, GC. Epidemiology and aetiology of heart failure. Nat. Rev. Cardiol.; 2016; 13, pp. 368-378. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26935038][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4868779][DOI: https://dx.doi.org/10.1038/nrcardio.2016.25]

3. Ahmad, F; Seidman, JG; Seidman, CE. The genetic basis for cardiac remodeling. Annu Rev. Genomics Hum. Genet; 2005; 6, pp. 185-216.1:CAS:528:DC%2BD2MXhtFWisrrL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16124859][DOI: https://dx.doi.org/10.1146/annurev.genom.6.080604.162132]

4. Dubois, CM et al. Evidence that furin is an authentic transforming growth factor-beta1-converting enzyme. Am. J. Pathol.; 2001; 158, pp. 305-316.1:CAS:528:DC%2BD3MXmsFOjtw%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11141505][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1850265][DOI: https://dx.doi.org/10.1016/S0002-9440(10)63970-3]

5. Nakajima, Y; Mironov, V; Yamagishi, T; Nakamura, H; Markwald, RR. Expression of smooth muscle alpha-actin in mesenchymal cells during formation of avian endocardial cushion tissue: a role for transforming growth factor beta3. Dev. Dyn.; 1997; 209, pp. 296-309.1:CAS:528:DyaK2sXkslCmtLo%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9215644][DOI: https://dx.doi.org/10.1002/(SICI)1097-0177(199707)209:3<296::AID-AJA5>3.0.CO;2-D]

6. Chen, MM; Lam, A; Abraham, JA; Schreiner, GF; Joly, AH. CTGF expression is induced by TGF- beta in cardiac fibroblasts and cardiac myocytes: a potential role in heart fibrosis. J. Mol. Cell Cardiol.; 2000; 32, pp. 1805-1819.1:CAS:528:DC%2BD3cXmvFWlu7k%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11013125][DOI: https://dx.doi.org/10.1006/jmcc.2000.1215]

7. Wang, B et al. The Kruppel-like factor KLF15 inhibits connective tissue growth factor (CTGF) expression in cardiac fibroblasts. J. Mol. Cell Cardiol.; 2008; 45, pp. 193-197.1:CAS:528:DC%2BD1cXhtVSjs7jI [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18586263][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2566509][DOI: https://dx.doi.org/10.1016/j.yjmcc.2008.05.005]

8. Campbell, SE; Katwa, LC. Angiotensin II stimulated expression of transforming growth factor-beta1 in cardiac fibroblasts and myofibroblasts. J. Mol. Cell Cardiol.; 1997; 29, pp. 1947-1958.1:CAS:528:DyaK2sXltVOqsrs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9236148][DOI: https://dx.doi.org/10.1006/jmcc.1997.0435]

9. Leask, A. Getting to the heart of the matter: new insights into cardiac fibrosis. Circ. Res; 2015; 116, pp. 1269-1276.1:CAS:528:DC%2BC2MXltlCktr4%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25814687][DOI: https://dx.doi.org/10.1161/CIRCRESAHA.116.305381]

10. Martin, P; Tzanidis, A; Stein-Oakley, A; Krum, H. Effect of a highly selective endothelin-converting enzyme inhibitor on cardiac remodeling in rats after myocardial infarction. J. Cardiovasc Pharm.; 2000; 36, pp. S367-S370.1:CAS:528:DC%2BD3cXnslKisb8%3D [DOI: https://dx.doi.org/10.1097/00005344-200036051-00106]

11. Porter, KE; Turner, NA. Cardiac fibroblasts: at the heart of myocardial remodeling. Pharm. Ther.; 2009; 123, pp. 255-278.1:CAS:528:DC%2BD1MXnsFyhtbY%3D [DOI: https://dx.doi.org/10.1016/j.pharmthera.2009.05.002]

12. Kong, P; Christia, P; Frangogiannis, NG. The pathogenesis of cardiac fibrosis. Cell Mol. Life Sci.; 2014; 71, pp. 549-574.1:CAS:528:DC%2BC2cXhtlSmsrc%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23649149][DOI: https://dx.doi.org/10.1007/s00018-013-1349-6]

13. Jackson, E. K. et al. SDF-1alpha (Stromal Cell-Derived Factor 1alpha) Induces cardiac fibroblasts, renal microvascular smooth muscle cells, and glomerular mesangial cells to proliferate, cause hypertrophy, and produce collagen. J Am Heart Assoc6, https://doi.org/10.1161/JAHA.117.007253 (2017).

14. Maurano, MT et al. Systematic localization of common disease-associated variation in regulatory DNA. Science; 2012; 337, pp. 1190-1195.2012Sci..337.1190M1:CAS:528:DC%2BC38Xht1ylsLfL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22955828][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3771521][DOI: https://dx.doi.org/10.1126/science.1222794]

15. Roadmap Epigenomics, C et al. Integrative analysis of 111 reference human epigenomes. Nature; 2015; 518, pp. 317-330. [DOI: https://dx.doi.org/10.1038/nature14248]

16. Umans, BD; Battle, A; Gilad, Y. Where are the disease-associated eQTLs?. Trends Genet; 2021; 37, pp. 109-124.1:CAS:528:DC%2BB3cXhslWqt7%2FL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32912663][DOI: https://dx.doi.org/10.1016/j.tig.2020.08.009]

17. Mostafavi, H., Spence, J. P., Naqvi, S. & Pritchard, J. K. Systematic differences in discovery of genetic effects on gene expression and complex traits. Nat. Genet55, 1866–1875 (2023).

18. Musunuru, K et al. From noncoding variant to phenotype via SORT1 at the 1p13 cholesterol locus. Nature; 2010; 466, pp. 714-719.2010Natur.466.714M1:CAS:528:DC%2BC3cXhtVSmurbO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20686566][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3062476][DOI: https://dx.doi.org/10.1038/nature09266]

19. Smemo, S et al. Obesity-associated variants within FTO form long-range functional connections with IRX3. Nature; 2014; 507, pp. 371-375.2014Natur.507.371S1:CAS:528:DC%2BC2cXks1Wlsbk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24646999][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4113484][DOI: https://dx.doi.org/10.1038/nature13138]

20. Gupta, RM et al. A genetic variant associated with five vascular diseases is a distal regulator of endothelin-1 gene expression. Cell; 2017; 170, pp. 522-533.e515.1:CAS:528:DC%2BC2sXht1GisLvL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28753427][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5785707][DOI: https://dx.doi.org/10.1016/j.cell.2017.06.049]

21. Wright, JB; Brown, SJ; Cole, MD. Upregulation of c-MYC in cis through a large chromatin loop linked to a cancer risk-associated single-nucleotide polymorphism in colorectal cancer cells. Mol. Cell Biol.; 2010; 30, pp. 1411-1420.1:CAS:528:DC%2BC3cXjsFaltrw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20065031][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2832500][DOI: https://dx.doi.org/10.1128/MCB.01384-09]

22. Cowper-Sal lari, R et al. Breast cancer risk-associated SNPs modulate the affinity of chromatin for FOXA1 and alter gene expression. Nat. Genet; 2012; 44, pp. 1191-1198.1:CAS:528:DC%2BC38Xhtlyls7jK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23001124][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3483423][DOI: https://dx.doi.org/10.1038/ng.2416]

23. Symmons, O et al. Functional and topological characteristics of mammalian regulatory domains. Genome Res; 2014; 24, pp. 390-400.1:CAS:528:DC%2BC2cXkvVyqs7s%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24398455][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3941104][DOI: https://dx.doi.org/10.1101/gr.163519.113]

24. Sexton, T; Cavalli, G. The role of chromosome domains in shaping the functional genome. Cell; 2015; 160, pp. 1049-1059.1:CAS:528:DC%2BC2MXks1agtbg%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25768903][DOI: https://dx.doi.org/10.1016/j.cell.2015.02.040]

25. Dixon, JR et al. Topological domains in mammalian genomes identified by analysis of chromatin interactions. Nature; 2012; 485, pp. 376-380.2012Natur.485.376D1:CAS:528:DC%2BC38XlslShurw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22495300][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3356448][DOI: https://dx.doi.org/10.1038/nature11082]

26. Schmitt, AD et al. A compendium of chromatin contact maps reveals spatially active regions in the human genome. Cell Rep.; 2016; 17, pp. 2042-2059.1:CAS:528:DC%2BC28XhvVOgsLjN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27851967][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5478386][DOI: https://dx.doi.org/10.1016/j.celrep.2016.10.061]

27. Rao, SS et al. A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping. Cell; 2014; 159, pp. 1665-1680.1:CAS:528:DC%2BC2cXitFCrtrrK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25497547][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5635824][DOI: https://dx.doi.org/10.1016/j.cell.2014.11.021]

28. Mifsud, B et al. Mapping long-range promoter contacts in human cells with high-resolution capture Hi-C. Nat. Genet; 2015; 47, pp. 598-606.1:CAS:528:DC%2BC2MXht1Wltb%2FL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25938943][DOI: https://dx.doi.org/10.1038/ng.3286]

29. Cairns, J et al. CHiCAGO: robust detection of DNA looping interactions in Capture Hi-C data. Genome Biol.; 2016; 17, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27306882][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4908757][DOI: https://dx.doi.org/10.1186/s13059-016-0992-2] 127.

30. Choy, MK et al. Promoter interactome of human embryonic stem cell-derived cardiomyocytes connects GWAS regions to cardiac gene networks. Nat. Commun.; 2018; 9, 2018NatCo..9.2526C [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29955040][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6023870][DOI: https://dx.doi.org/10.1038/s41467-018-04931-0] 2526.

31. Montefiori, L. E. et al. A promoter interaction map for cardiovascular disease genetics. Elifehttps://doi.org/10.7554/eLife.35788 (2018).

32. Arvanitis, M et al. Genome-wide association and multi-omic analyses reveal ACTN2 as a gene linked to heart failure. Nat. Commun.; 2020; 11, 2020NatCo.11.1122A [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32111823][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7048760][DOI: https://dx.doi.org/10.1038/s41467-020-14843-7] 1122.

33. Anene-Nzelu, CG et al. Assigning distal genomic enhancers to cardiac disease-causing genes. Circulation; 2020; 142, pp. 910-912. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32866060][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9136986][DOI: https://dx.doi.org/10.1161/CIRCULATIONAHA.120.046040]

34. Hocker, J. D. et al. Cardiac cell type-specific gene regulatory programs and disease risk association. Sci. Adv.https://doi.org/10.1126/sciadv.abf1444 (2021).

35. Andersson, R et al. An atlas of active enhancers across human cell types and tissues. Nature; 2014; 507, pp. 455-461.2014Natur.507.455A1:CAS:528:DC%2BC2cXkvV2mt7g%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24670763][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5215096][DOI: https://dx.doi.org/10.1038/nature12787]

36. Christophersen, IE et al. Large-scale analyses of common and rare variants identify 12 new loci associated with atrial fibrillation. Nat. Genet; 2017; 49, pp. 946-952.1:CAS:528:DC%2BC2sXmtFGhtL0%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28416818][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5585859][DOI: https://dx.doi.org/10.1038/ng.3843]

37. Nielsen, JB et al. Biobank-driven genomic discovery yields new insight into atrial fibrillation biology. Nat. Genet; 2018; 50, pp. 1234-1239.1:CAS:528:DC%2BC1cXhsVWru7rF [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30061737][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6530775][DOI: https://dx.doi.org/10.1038/s41588-018-0171-3]

38. van der Harst, P; Verweij, N. Identification of 64 Novel Genetic Loci Provides an Expanded View on the Genetic Architecture of Coronary Artery Disease. Circ. Res; 2018; 122, pp. 433-443. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29212778][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5805277][DOI: https://dx.doi.org/10.1161/CIRCRESAHA.117.312086]

39. Zhang, S; Ding, Y; Feng, F; Gao, Y. The role of blood CXCL12 level in prognosis of coronary artery disease: A meta-analysis. Front Cardiovasc Med; 2022; 9, 1:CAS:528:DC%2BB38XisVWjs7zM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35966557][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9363627][DOI: https://dx.doi.org/10.3389/fcvm.2022.938540] 938540.

40. Nikpay, M et al. A comprehensive 1000 Genomes-based genome-wide association meta-analysis of coronary artery disease. Nat. Genet; 2015; 47, pp. 1121-1130.1:CAS:528:DC%2BC2MXhsV2lt7zO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26343387][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4589895][DOI: https://dx.doi.org/10.1038/ng.3396]

41. Nelson, CP et al. Association analyses based on false discovery rate implicate new loci for coronary artery disease. Nat. Genet; 2017; 49, pp. 1385-1391.1:CAS:528:DC%2BC2sXhtFOru7vE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28714975][DOI: https://dx.doi.org/10.1038/ng.3913]

42. Amrute, JM et al. Targeting immune-fibroblast cell communication in heart failure. Nature; 2024; 635, pp. 423-433.2024Natur.635.423A1:CAS:528:DC%2BB2cXitlanurrO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39443792][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC12334188][DOI: https://dx.doi.org/10.1038/s41586-024-08008-5]

43. Kuppe, C et al. Spatial multi-omic map of human myocardial infarction. Nature; 2022; 608, pp. 766-777.2022Natur.608.766K1:CAS:528:DC%2BB38XitFSjtb3O [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35948637][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9364862][DOI: https://dx.doi.org/10.1038/s41586-022-05060-x]

44. Amrute, J. M. et al. Targeting Immune-Fibroblast Crosstalk in Myocardial Infarction and Cardiac Fibrosis. Res Sq, https://doi.org/10.21203/rs.3.rs-2402606/v1 (2023).

45. Abugessaisa, I; Noguchi, S; Carninci, P; Kasukawa, T. The FANTOM5 Computation Ecosystem: Genomic Information Hub for Promoters and Active Enhancers. Methods Mol. Biol.; 2017; 1611, pp. 199-217.1:CAS:528:DC%2BC1cXhsFWnt7vO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28451981][DOI: https://dx.doi.org/10.1007/978-1-4939-7015-5_15]

46. Forcato, M et al. Comparison of computational methods for Hi-C data analysis. Nat. Methods; 2017; 14, pp. 679-685.1:CAS:528:DC%2BC2sXpvFCmsL4%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28604721][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5493985][DOI: https://dx.doi.org/10.1038/nmeth.4325]

47. Zhu, Z et al. Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets. Nat. Genet; 2016; 48, pp. 481-487.1:CAS:528:DC%2BC28XkvFGrs7g%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27019110][DOI: https://dx.doi.org/10.1038/ng.3538]

48. Miyazawa, K; Ito, K. The Evolving Story in the Genetic Analysis for Heart Failure. Front Cardiovasc Med; 2021; 8, 1:CAS:528:DC%2BB3MXitFClsbzJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33928132][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8076510][DOI: https://dx.doi.org/10.3389/fcvm.2021.646816] 646816.

49. Asazuma-Nakamura, Y et al. Cx43 contributes to TGF-beta signaling to regulate differentiation of cardiac fibroblasts into myofibroblasts. Exp. Cell Res; 2009; 315, pp. 1190-1199.1:CAS:528:DC%2BD1MXktF2ksbY%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19162006][DOI: https://dx.doi.org/10.1016/j.yexcr.2008.12.021]

50. Panek, AN et al. Connective tissue growth factor overexpression in cardiomyocytes promotes cardiac hypertrophy and protection against pressure overload. PLoS One; 2009; 4, e6743.2009PLoSO..4.6743P [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19707545][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2727794][DOI: https://dx.doi.org/10.1371/journal.pone.0006743]

51. Li, Y et al. Global genetic analysis in mice unveils central role for cilia in congenital heart disease. Nature; 2015; 521, pp. 520-524.2015Natur.521.520L1:CAS:528:DC%2BC2MXlsFektLk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25807483][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4617540][DOI: https://dx.doi.org/10.1038/nature14269]

52. Huttemann, M et al. Mice deleted for heart-type cytochrome c oxidase subunit 7a1 develop dilated cardiomyopathy. Mitochondrion; 2012; 12, pp. 294-304.1:CAS:528:DC%2BC38XjsFKrtLw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22119795][DOI: https://dx.doi.org/10.1016/j.mito.2011.11.002]

53. Doring, Y; Pawig, L; Weber, C; Noels, H. The CXCL12/CXCR4 chemokine ligand/receptor axis in cardiovascular disease. Front Physiol.; 2014; 5, 212. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24966838][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4052746]

54. Li, Y et al. Cardiac Fibroblast-Specific Activating Transcription Factor 3 Protects Against Heart Failure by Suppressing MAP2K3-p38 Signaling. Circulation; 2017; 135, pp. 2041-2057.1:CAS:528:DC%2BC2sXot1yntLs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28249877][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5542579][DOI: https://dx.doi.org/10.1161/CIRCULATIONAHA.116.024599]

55. Shi, Y; Riese, DJ, 2nd; Shen, J. The Role of the CXCL12/CXCR4/CXCR7 Chemokine Axis in Cancer. Front Pharm.; 2020; 11, 574667.1:CAS:528:DC%2BB3MXjslyrur8%3D [DOI: https://dx.doi.org/10.3389/fphar.2020.574667]

56. Tanaka, Y et al. Systems analysis of ATF3 in stress response and cancer reveals opposing effects on pro-apoptotic genes in p53 pathway. PLoS One; 2011; 6, e26848.2011PLoSO..626848T1:CAS:528:DC%2BC3MXhsVKlur7N [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22046379][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3202577][DOI: https://dx.doi.org/10.1371/journal.pone.0026848]

57. Collaboration, I. R. G. C. E. R. F.et al. Interleukin-6 receptor pathways in coronary heart disease: a collaborative meta-analysis of 82 studies. Lancet; 2012; 379, pp. 1205-1213. [DOI: https://dx.doi.org/10.1016/S0140-6736(11)61931-4]

58. Gonzalez, GE et al. Deletion of interleukin-6 prevents cardiac inflammation, fibrosis and dysfunction without affecting blood pressure in angiotensin II-high salt-induced hypertension. J. Hypertens.; 2015; 33, pp. 144-152.1:CAS:528:DC%2BC2cXitVamtL%2FJ [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25304471][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4256143][DOI: https://dx.doi.org/10.1097/HJH.0000000000000358]

59. Interleukin-6 Receptor Mendelian Randomisation Analysis, C.et al. The interleukin-6 receptor as a target for prevention of coronary heart disease: a mendelian randomisation analysis. Lancet; 2012; 379, pp. 1214-1224. [DOI: https://dx.doi.org/10.1016/S0140-6736(12)60110-X]

60. Saadane, N; Alpert, L; Chalifour, LE. Altered molecular response to adrenoreceptor-induced cardiac hypertrophy in Egr-1-deficient mice. Am. J. Physiol. Heart Circ. Physiol.; 2000; 278, pp. H796-H805.1:CAS:528:DC%2BD3cXhvFaqur8%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/10710348][DOI: https://dx.doi.org/10.1152/ajpheart.2000.278.3.H796]

61. Liu, Y et al. Snail1 is involved in de novo cardiac fibrosis after myocardial infarction in mice. Acta Biochim Biophys. Sin. (Shanghai); 2012; 44, pp. 902-910.1:CAS:528:DC%2BC38Xhs1Sqs7%2FI [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23059020][DOI: https://dx.doi.org/10.1093/abbs/gms085]

62. Gamazon, ER et al. Using an atlas of gene regulation across 44 human tissues to inform complex disease- and trait-associated variation. Nat. Genet; 2018; 50, pp. 956-967.1:CAS:528:DC%2BC1cXht1Cjt7jK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29955180][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6248311][DOI: https://dx.doi.org/10.1038/s41588-018-0154-4]

63. Ichiki, T et al. Differential expression of the pro-natriuretic peptide convertases corin and furin in experimental heart failure and atrial fibrosis. Am. J. Physiol. Regul. Integr. Comp. Physiol.; 2013; 304, pp. R102-R109.1:CAS:528:DC%2BC3sXit1ygt7s%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23152112][DOI: https://dx.doi.org/10.1152/ajpregu.00233.2012]

64. Yakala, GK et al. FURIN Inhibition Reduces Vascular Remodeling and Atherosclerotic Lesion Progression in Mice. Arterioscler Thromb. Vasc. Biol.; 2019; 39, pp. 387-401.1:CAS:528:DC%2BC1MXjslKktLs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30651003][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6393193][DOI: https://dx.doi.org/10.1161/ATVBAHA.118.311903]

65. Vandekerckhove, J; Bugaisky, G; Buckingham, M. Simultaneous expression of skeletal muscle and heart actin proteins in various striated muscle tissues and cells. A quantitative determination of the two actin isoforms. J. Biol. Chem.; 1986; 261, pp. 1838-1843.1:CAS:528:DyaL28XhtVKrtr8%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/3944112][DOI: https://dx.doi.org/10.1016/S0021-9258(17)36017-9]

66. Fairfax, BP et al. Innate immune activity conditions the effect of regulatory variants upon monocyte gene expression. Science; 2014; 343, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24604202][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4064786][DOI: https://dx.doi.org/10.1126/science.1246949] 1246949.

67. Farh, KK et al. Genetic and epigenetic fine mapping of causal autoimmune disease variants. Nature; 2015; 518, pp. 337-343.2015Natur.518.337F1:CAS:528:DC%2BC2cXhvVemtrbP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25363779][DOI: https://dx.doi.org/10.1038/nature13835]

68. Replogle, JM et al. Mapping information-rich genotype-phenotype landscapes with genome-scale Perturb-seq. Cell; 2022; 185, pp. 2559-2575.e2528.1:CAS:528:DC%2BB38XhsFWgsLrP [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35688146][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9380471][DOI: https://dx.doi.org/10.1016/j.cell.2022.05.013]

69. Metzner, E., Southard, K. M. & Norman, T. M. Multiome Perturb-seq unlocks scalable discovery of integrated perturbation effects on the transcriptome and epigenome. Cell Syst16, 101161, https://doi.org/10.1016/j.cels.2024.12.002 (2025).

70. Servant, N et al. HiC-Pro: an optimized and flexible pipeline for Hi-C data processing. Genome Biol.; 2015; 16, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26619908][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4665391][DOI: https://dx.doi.org/10.1186/s13059-015-0831-x] 259.

71. Durand, NC et al. Juicer Provides a One-Click System for Analyzing Loop-Resolution Hi-C Experiments. Cell Syst.; 2016; 3, pp. 95-98.1:CAS:528:DC%2BC2sXhtFKksbk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27467249][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5846465][DOI: https://dx.doi.org/10.1016/j.cels.2016.07.002]

72. Mifsud, B et al. GOTHiC, a probabilistic model to resolve complex biases and to identify real interactions in Hi-C data. PLoS One; 2017; 12, e0174744. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28379994][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5381888][DOI: https://dx.doi.org/10.1371/journal.pone.0174744]

73. Schoenfelder, S et al. The pluripotent regulatory circuitry connecting promoters to their long-range interacting elements. Genome Res; 2015; 25, pp. 582-597.1:CAS:528:DC%2BC2MXmvFCgu7c%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25752748][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4381529][DOI: https://dx.doi.org/10.1101/gr.185272.114]

74. Lee, D-I; Roy, S. GRiNCH: simultaneous smoothing and detection of topological units of genome organization from sparse chromatin contact count matrices with matrix factorization. Genome Biol.; 2021; 22, 1:CAS:528:DC%2BB3MXht1KrtLzI [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34034791][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8152090][DOI: https://dx.doi.org/10.1186/s13059-021-02378-z] 164.

75. Kruse, K; Hug, CB; Hernández-Rodríguez, B; Vaquerizas, JM. TADtool: visual parameter identification for TAD-calling algorithms. Bioinformatics; 2016; 32, pp. 3190-3192.1:CAS:528:DC%2BC2sXht1ajtbrK [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27318199][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5048066][DOI: https://dx.doi.org/10.1093/bioinformatics/btw368]

76. Quinlan, AR; Hall, IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics; 2010; 26, pp. 841-842.1:CAS:528:DC%2BC3cXivFGkurc%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20110278][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2832824][DOI: https://dx.doi.org/10.1093/bioinformatics/btq033]

77. Buenrostro, JD; Wu, B; Chang, HY; Greenleaf, WJ. ATAC-seq: A Method for Assaying Chromatin Accessibility Genome-Wide. Curr. Protoc. Mol. Biol.; 2015; 109, pp. 21 29 21-21 29 29. [DOI: https://dx.doi.org/10.1002/0471142727.mb2129s109]

78. Lee, J. W., Chuan, S. F., Kim, D., Boley, N., Kundaje, A. ATAC-Seq / DNase-Seq Pipeline, https://github.com/kundajelab/atac_dnase_pipelines (2017).

79. Gusmao, EG; Allhoff, M; Zenke, M; Costa, IG. Analysis of computational footprinting methods for DNase sequencing experiments. Nat. Methods; 2016; 13, pp. 303-309. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26901649][DOI: https://dx.doi.org/10.1038/nmeth.3772]

80. Gusmao, EG; Dieterich, C; Zenke, M; Costa, IG. Detection of active transcription factor binding sites with the combination of DNase hypersensitivity and histone modifications. Bioinformatics; 2014; 30, pp. 3143-3151.1:CAS:528:DC%2BC28XhtFOrtLnN [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25086003][DOI: https://dx.doi.org/10.1093/bioinformatics/btu519]

81. Li, Z et al. Identification of transcription factor binding sites using ATAC-seq. Genome Biol.; 2019; 20, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30808370][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6391789][DOI: https://dx.doi.org/10.1186/s13059-019-1642-2] 45.

82. Li, Z et al. RGT: a toolbox for the integrative analysis of high throughput regulatory genomics data. BMC Bioinforma.; 2023; 24, 1:CAS:528:DC%2BB3sXlsV2qtr0%3D [DOI: https://dx.doi.org/10.1186/s12859-023-05184-5] 79.

83. Kulakovskiy, IV et al. HOCOMOCO: towards a complete collection of transcription factor binding models for human and mouse via large-scale ChIP-Seq analysis. Nucleic Acids Res; 2018; 46, pp. D252-D259.1:CAS:528:DC%2BC1cXitlGitr%2FE [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29140464][DOI: https://dx.doi.org/10.1093/nar/gkx1106]

84. Khan, A et al. JASPAR 2018: update of the open-access database of transcription factor binding profiles and its web framework. Nucleic Acids Res; 2018; 46, pp. D260-D266.1:CAS:528:DC%2BC1cXitlGitr3N [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29140473][DOI: https://dx.doi.org/10.1093/nar/gkx1126]

85. Sandelin, A; Alkema, W; Engstrom, P; Wasserman, WW; Lenhard, B. JASPAR: an open-access database for eukaryotic transcription factor binding profiles. Nucleic Acids Res; 2004; 32, pp. D91-D94.1:CAS:528:DC%2BD3sXhtVSrur7N [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/14681366][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC308747][DOI: https://dx.doi.org/10.1093/nar/gkh012]

86. Wilczynski, B; Dojer, N; Patelak, M; Tiuryn, J. Finding evolutionarily conserved cis-regulatory modules with a universal set of motifs. BMC Bioinforma.; 2009; 10, [DOI: https://dx.doi.org/10.1186/1471-2105-10-82] 82.

87. Heinz, S et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell; 2010; 38, pp. 576-589.1:CAS:528:DC%2BC3cXns1SlsLc%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20513432][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2898526][DOI: https://dx.doi.org/10.1016/j.molcel.2010.05.004]

88. Parkhomchuk, D et al. Transcriptome analysis by strand-specific sequencing of complementary DNA. Nucleic Acids Res; 2009; 37, e123. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19620212][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2764448][DOI: https://dx.doi.org/10.1093/nar/gkp596]

89. Sultan, M et al. A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome. Science; 2008; 321, pp. 956-960.2008Sci..321.956S1:CAS:528:DC%2BD1cXpslWrur4%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18599741][DOI: https://dx.doi.org/10.1126/science.1160342]

90. MacArthur, J et al. The new NHGRI-EBI Catalog of published genome-wide association studies (GWAS Catalog. Nucleic Acids Res; 2017; 45, pp. D896-D901.1:CAS:528:DC%2BC1cXhslWhs7c%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27899670][DOI: https://dx.doi.org/10.1093/nar/gkw1133]

91. Shah, S et al. Genome-wide association and Mendelian randomisation analysis provide insights into the pathogenesis of heart failure. Nat. Commun.; 2020; 11, 2020NatCo.11.163S1:CAS:528:DC%2BB3cXpt1yntg%3D%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31919418][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6952380][DOI: https://dx.doi.org/10.1038/s41467-019-13690-5] 163.

92. Machiela, MJ; Chanock, SJ. LDlink: a web-based application for exploring population-specific haplotype structure and linking correlated alleles of possible functional variants. Bioinformatics; 2015; 31, pp. 3555-3557.1:CAS:528:DC%2BC28XhsVersbfL [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26139635][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4626747][DOI: https://dx.doi.org/10.1093/bioinformatics/btv402]

93. Yates, B et al. Genenames.org: the HGNC and VGNC resources in 2017. Nucleic Acids Res; 2017; 45, pp. D619-D625.1:CAS:528:DC%2BC1cXhslWhtLg%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27799471][DOI: https://dx.doi.org/10.1093/nar/gkw1033]

94. Carithers, LJ et al. A novel approach to high-quality postmortem tissue procurement: the GTEx project. Biopreserv Biobank; 2015; 13, pp. 311-319. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26484571][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4675181][DOI: https://dx.doi.org/10.1089/bio.2015.0032]

95. Durand, NC et al. Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom. Cell Syst.; 2016; 3, pp. 99-101.1:CAS:528:DC%2BC2sXhtFKks7w%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27467250][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5596920][DOI: https://dx.doi.org/10.1016/j.cels.2015.07.012]

96. Kerpedjiev P. et al. HiGlass: web-based visual exploration and analysis of genome interaction maps. Genome Biol.19, 125 (2018).

97. Kramer, NE et al. Plotgardener: cultivating precise multi-panel figures in R. Bioinformatics; 2022; 38, pp. 2042-2045.1:CAS:528:DC%2BB38XhtFygs7jM [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35134826][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8963281][DOI: https://dx.doi.org/10.1093/bioinformatics/btac057]

98. Sigma, M. CRISPR Design Tools, https://www.milliporesigmabioinfo.com/bioinfo_tools/faces/secured/crispr/crispr.xhtml (2021).

99. Labun, K et al. CHOPCHOP v3: expanding the CRISPR web toolbox beyond genome editing. Nucleic Acids Res; 2019; 47, pp. W171-W174.1:CAS:528:DC%2BB3cXktVyiu7c%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31106371][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6602426][DOI: https://dx.doi.org/10.1093/nar/gkz365]

100. Pliatsika, V; Rigoutsos, I. “Off-Spotter”: very fast and exhaustive enumeration of genomic lookalikes for designing CRISPR/Cas guide RNAs. Biol. Direct; 2015; 10, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25630343][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4326336][DOI: https://dx.doi.org/10.1186/s13062-015-0035-z] 4.

101. Bae, S; Park, J; Kim, JS. Cas-OFFinder: a fast and versatile algorithm that searches for potential off-target sites of Cas9 RNA-guided endonucleases. Bioinformatics; 2014; 30, pp. 1473-1475.1:CAS:528:DC%2BC2cXnslOitrs%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24463181][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4016707][DOI: https://dx.doi.org/10.1093/bioinformatics/btu048]

102. Hwang, GH; Kim, JS; Bae, S. Web-Based CRISPR Toolkits: Cas-OFFinder, Cas-Designer, and Cas-Analyzer. Methods Mol. Biol.; 2021; 2162, pp. 23-33.1:CAS:528:DC%2BB3MXptVems7Y%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32926375][DOI: https://dx.doi.org/10.1007/978-1-0716-0687-2_2]

103. Horlbeck, M. A. et al. Compact and highly active next-generation libraries for CRISPR-mediated gene repression and activation. Elifehttps://doi.org/10.7554/eLife.19760 (2016).

104. Satpathy, AT et al. Massively parallel single-cell chromatin landscapes of human immune cell development and intratumoral T cell exhaustion. Nat. Biotechnol.; 2019; 37, pp. 925-936.1:CAS:528:DC%2BC1MXhsFWqt7fO [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31375813][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7299161][DOI: https://dx.doi.org/10.1038/s41587-019-0206-z]

105. Mimitou, EP et al. Multiplexed detection of proteins, transcriptomes, clonotypes and CRISPR perturbations in single cells. Nat. Methods; 2019; 16, pp. 409-412.1:CAS:528:DC%2BC1MXosFKmurk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31011186][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6557128][DOI: https://dx.doi.org/10.1038/s41592-019-0392-0]

106. Stuart, T et al. Comprehensive Integration of single-cell data. Cell; 2019; 177, pp. 1888-1902.e1821.1:CAS:528:DC%2BC1MXhtFens77L [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31178118][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6687398][DOI: https://dx.doi.org/10.1016/j.cell.2019.05.031]

107. Tirosh, I et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science; 2016; 352, pp. 189-196.2016Sci..352.189T1:CAS:528:DC%2BC28XlsFSmurk%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27124452][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4944528][DOI: https://dx.doi.org/10.1126/science.aad0501]

108. Papalexi, E et al. Characterizing the molecular regulation of inhibitory immune checkpoints with multimodal single-cell screens. Nat. Genet; 2021; 53, pp. 322-331.1:CAS:528:DC%2BB3MXlsV2lt7Y%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33649593][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8011839][DOI: https://dx.doi.org/10.1038/s41588-021-00778-2]

109. Finak, G et al. MAST: a flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol.; 2015; 16, [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26653891][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4676162][DOI: https://dx.doi.org/10.1186/s13059-015-0844-5] 278.

110. Wu, T et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innov. (Camb.); 2021; 2, 100141.1:CAS:528:DC%2BB38XitVemt7jI

111. Yu, G; Wang, LG; Han, Y; He, QY. ClusterProfiler: an R package for comparing biological themes among gene clusters. OMICS; 2012; 16, pp. 284-287.1:CAS:528:DC%2BC38XmsFarsLw%3D [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22455463][PubMedCentral: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3339379][DOI: https://dx.doi.org/10.1089/omi.2011.0118]