Patients and controls included in the exome-wide association analysis have been reported in our previous study (Xie et al., 2019). Briefly, the complete cohort of patients with complex VSD comprised 60 PA/VSD cases and 20 TOF cases, whereas the control cohort consisted of 100 healthy non-hospitalized children. Written informed consents were obtained from all participants or their guardians before enrollment. This study compiled with the Helsinki Declaration and was approved by the Ethics Committee of Xinhua Hospital affiliated to Shanghai Jiao Tong University (No. XHEC-C-2012-018).

WES data of 80 complex VSD cases and 100 controls were derived from our previous study (Xie et al., 2019). No individuals were excluded because of low genotyping rate (<95%) or sex discrepancies. SNPs were excluded if they: 1) were not located on autosomal chromosomes; 2) had a call rate of <95%; 3) had a MAF of <0.05 in controls; 4) had a p value of <1 × 10⁻⁵ in the Hardy-Weinberg equilibrium test among controls. A set of 125,399 SNPs passed quality control criteria and was used for exome-wide association analysis.

Population structure was evaluated by principal component analysis (PCA) as implemented via PLINK (version 1.90) (Purcell et al., 2007). The first two eigenvectors were visualized via R package “ggplot2”. Exome-wide associations were assessed in an additive model by logistic regression analysis with adjustment for the top eigenvector. SNPs that had a p value below 5 × 10⁻⁸ were considered exome-wide significant.

Results of exome-wide association analysis were functional annotated via FUMA (http://fuma.ctglab.nl/) (Watanabe et al., 2017). Exome-wide significant SNPs that were independent from each other at r ² < 0.6 were defined as independent lead SNPs. All filtered SNPs that were in linkage disequilibrium (LD) (r ² > 0.2) with one of the independent lead SNPs were included in the LD expansion. Subsequently, causal SNPs were defined by merging independent lead SNPs and corresponding LD SNPs. Functional annotation for these SNPs were acquired via ANNOVAR using Ensembl genes (Wang K. et al., 2010).

ATAC-seq libraries in human embryonic stem cells (hESCs) were generated by Qing et al. (2019), and in hESC-derived cardiac lineage cells generated by Bertero et al. (2019). Raw fastq files were downloaded from GSE109524 for hESCs, and from GSE106690 for cardiac lineage cells. ATAC-seq reads were trimmed via Trim_galore (version 3.4) and mapped to the hg19 genome via Bowtie2 (version 2.4.4) with default parameters. Uniquely mapping reads were kept by filtering for their primary mapping location with MAPQ score >25. Duplicate reads were removed via Samtools (version 1.9). Significant peaks were called via MACS2 (version 2.1.1.20160309) using the parameters “--nomodel--shift-100 --extsize 200 -q 0.05”. Gene tracks were normalized to reads per genomic content (RPGC) and visualized via Integrative Genomics Viewers (IGV).

Genome-wide correlation heatmap was generated using pearson correlation method via DeepTools (version 3.4.3) multiBamSummary (Ramirez et al., 2016). Processed bam files were used as inputs. All chromosomes were divided into 10 kb non-overlapping bins in the correlation analysis. The resulting matrix was visualized via R package “pheatmap”.

Read density values were calculated via DeepTools computeMatrix using the parameters “--referencePoint TSS -a 1000 -b 1000”. The resulting matrices were visualized via DeepTools plotHeatmap. Of note, read density heatmaps were sorted by mean coverage per locus.

H3K27ac ChIP-seq libraries of normal human heart tissues [specifically, right atrium (RA), left ventricle (LV) and right ventricle (RV)] were generated by the NIH Roadmap Epigenomics Mapping Consortium (Bernstein et al., 2010). Raw fastq files were downloaded from GSE16256. Analysis pipeline for ChIP-seq reads was the same as that for ATAC-seq. Significant narrow peaks were called via MACS2 with default parameters. Gene tracks were normalized to RPGC and visualized via IGV.

In situ promoter CHi-C libraries in human induced pluripotent stem cells (iPSCs) and iPSC-derived CMs were generated by Montefiori et al. (2018). Processed sequencing data were downloaded from E-MTAB-6014 to find interactions between SNPs and promoters. Selected interactions were visualized via IGV.

Microarray data of human embryonic heart samples from Carnegie stages 11–15 were obtained from our former study (Shi et al., 2018). Gene expression data of human adult hearts were retrieved from the Genotype-Tissue Expression (GTEx) database (GTEx Consortium, 2020).

We proposed an integrated multi-omic analysis to systematically recognize triplets of cis-regulatory SNP, relevant TF and target gene in complex VSD. Our hypothesis was that cis-regulatory SNPs modulating the risk of complex VSD would probably reside in cis-regulatory regions to disrupt the binding affinity of corresponding TFs, which then could alter transcription activity of their target genes in the heart. The schematic illustration of our framework was shown in Figure 1. Briefly, “cis-regulatory SNP-TF” pairs were identified by combining exome-wide association analysis, epigenetic annotation of cis-regulatory regions and motif analysis. Putative target genes were mainly determined based on promoter CHi-C analysis and binding site prediction.

Population structure of the 80 complex VSD cases and 100 controls was described in PCA plot (Figure 2B). There was remarkable population stratification between cases and controls, but little heterogeneity between PA/VSD and TOF individuals. To reduce false associations induced by population stratification, we assessed the exome-wide association in an additive logistic regression model with adjustment for the top eigenvector. The overall p values were presented in Figure 2A and 266 SNPs exhibited a significant association with p < 5 × 10⁻⁸ (Supplementary Table S1). Of all the exome-wide significant SNPs, 95 (35%) were in exonic regions, 133 (50%) were in intronic regions, 12 (5%) were in UTR5/UTR3 regions, and 26 (10%) were in intergenic regions (Figure 2C). We used FUMA to functionally annotate results from the association analysis and extracted 91 independent lead SNPs. Noteworthy, most independent lead SNPs in the locus were not necessarily the causal SNP but that it might be in LD with the causal one (Bush and Moore, 2012). Following the criteria of r ² > 0.2, we examined every 2-Mb region centered on each lead SNP and identified 143 LD SNPs. Together, our exome-wide association analysis obtained a total of 194 causal SNPs (Supplementary Table S2). Of all the causal SNPs, 140 (72%) were in non-coding regions, and thereby were included in subsequent functional analysis for cis-regulatory SNP characterization (Figure 2C).

Cardiac morphogenesis is a tightly-regulated process involving multiple stages and perturbations in any stage can lead to complex VSD. In vitro directed differentiation of hESC into CM serves as an ideal system to elucidate spatiotemporal regulation of cardiogenesis (Murry and Keller, 2008). As chromatin accessibility was a hallmark of genomic regulatory regions, we analyzed public ATAC-seq datasets for hESC, mesodermal cell (MES), cardiac progenitor cell (CPC) and CM to obtain cardiac-specific cis-regulatory regions (Figure 3A). Interestingly, chromatin accessibility patterns of MES, CPC, and CM were extremely similar as revealed by ATAC-seq, but these patterns were far different from that of hESC (Figures 3B,C), indicating a strong conservation of cis-regulatory regions within cardiac lineage cells. Considering the previous report that such regulatory regions were cell-type specific in brain (Nott et al., 2019), we next asked whether our finding reflected the nature of cardiac morphogenesis. A high-resolution examination of ATAC-seq enrichment near lineage specification genes was conducted in all in vitro cardiogenic stages, which suggested a high correlation between chromatin accessibility and transcriptional activity. Notable examples included NANOG in the hESC, PDGFRA in the MES, GATA4 in the CPC and MYH7 in the CM (Figures 3D–G). Significant ATAC-seq peaks were observed at promoters or putative enhancers at the same stages when they were expressed. Thus, owing to the similar cardiogenic expression pattern in distinct cardiac lineage cells, it is convincing that chromatin accessibility and cis-regulatory regions remain little changed during cardiac differentiation process.

According to the retrieved ATAC-seq datasets, there were 99,666 open chromatin regions in MES, 82,810 in CPC as well as 48,278 in CM (Figure 3H). A total of 134,843 open chromatin regions were identified during cardiac morphogenesis after merging those in MES, CPC, and CM, and were recognized as cardiac-specific cis-regulatory regions.

The schematic representation of identifying cis-regulatory SNP candidates for complex VSD was shown in Figure 3I. We found three of 140 causal non-coding SNPs resided in cardiac-specific cis-regulatory regions, namely rs2279658, rs4794004, and rs2291726 (Figure 3I). Of note, the SNP rs2279658 was located at the center of a short sequence conserved among humans, rhesus monkeys and elephants (Figure 4A). On the contrary, genomic sequences containing rs4794004-A or rs2291726-T varied among the aforementioned species (Figures 4B,C). Given that these three SNPs lay in intergenic regions, we then asked whether their corresponding cis-regulatory regions had enhancer activity in the heart. Histone acetylation such as H3K27ac has been well recognized as a marker for active enhancers and an excellent indicator of enhancer activity (Creyghton et al., 2010). Thus, we retrieved H3K27ac ChIP-seq datasets produced by the NIH Roadmap Epigenomics Mapping Consortium from normal human heart tissues of left atrium as well as left and right ventricle (Bernstein et al., 2010). There was significant H3K27ac enrichment at the rs2279658-A locus in all tissues examined, whereas the rest two exhibited undetectable H3K27ac signals (Figures 4D,E). With these above evidences we inferred rs2279658 to be an important cis-regulatory candidate associated with complex VSD.

To predict TFs potentially affected by cis-regulatory SNP rs2279658, we profiled direct interactions between the rs2279658-A containing sequence and TFs via JASPAR database (Castro-Mondragon et al., 2022). Several well-studied TFs were found to recognize and bind to the selected sequence, including CRX, FOXH1, GSC, GSC2, OTX2, PITX1, and PITX2 (Figure 4F). Among these TF candidates, FOXH1 and PITX2 have been formerly reported to play critical roles in embryonic heart development and thereby triggered our interest. Motif analysis revealed that both FOXH1 and PITX2 recognized a 5-bp consensus sequence “AATCC” to promote or inhibit the assembly of transcription machinery (Figure 4G), so that the risk allele rs2279658-G would abrogate binding of these 2 TFs to corresponding cis-regulatory regions. Besides, considerable mRNA levels of FOXH1 and PITX2 were observed in human embryonic heart from Carnegie stage 11–15 (Figure 4H), despite few in human adult heart (Supplementary Figures S1A,B). In general, the risk allele rs2279658-G can indeed be associated with complex VSD in a cis-regulatory manner.

We explored a public promoter interaction data for cardiovascular disease genetics to pinpoint target genes regulated by rs2279658. The chromatin interaction map, developed via in situ promoter CHi-C recorded high-resolution promoter-enhancer interactions of human iPSCs and iPSC-derived CMs (Montefiori et al., 2018). rs2279658 was found to link to the PLAC8 gene promoter in iPSCs (p < 0.05) (Figure 5A), and to LIN54, COQ2, and FAM175A gene in CMs (p < 0.05) (Figure 5B). This data accorded with the above-mentioned notion that cis-regulatory regions varied among cell types, so that such promoter-enhancer interactome was applicable for profiling “cis-regulatory SNP-target gene” pairs in complex VSD. Since the CM promoter interaction map linked SNPs to target genes relevant with cardiovascular disease, LIN54, COQ2, and FAM175A were prioritized for further TF correlation analysis.

According to the JASPAR database, there existed two binding sites of FOXH1 in the COQ2 gene promoter (Figure 5C), as well as binding sites of both FOXH1 and PITX2 in the FAM175A promoter (Figure 5C). By contrast, the binding site of neither FOXH1 nor PITX2 was predicted in the LIN54 gene promoter (data not shown). Tissue expression analysis showed that expression of COQ2 and FAM175A could be detected in fetal heart from Carnegie stage 11–15 (Figure 5D) as well as human adult heart (Supplementary Figures S1C,D). The risk allele rs2279658-G was also related to decreased expression of FAM175A in the whole blood (data not shown) (Jansen et al., 2017). Taken together, we conjectured that the rs2279658 A>G change might affect complex VSD risk by changing the level of COQ2 or FAM175A.