Bioinformatics and In-Silico Findings Reveal Candidate Genes for Tetralogy of Fallot via Integrative Multi-Omics Data

1 Introduction

Tetralogy of Fallot (TOF), the most common cyanotic congenital heart disease, is characterized by four cardinal anatomical abnormalities: ventricular septal defects, pulmonary artery stenosis, aortic cross, and right ventricular hypertrophy [1,2]. Epidemiological data have shown that TOF occurs in 3–5 per 10,000 newborns and the family recurrence risk of TOF approximating 3% [3,4,5]. Advances in congenital cardiac surgical techniques have significantly reduced early mortality rates from 25% to less than 2% in recent decades [6]. However, postoperative morbidity and mortality remain substantial for patients with complex congenital heart disease, particularly among infants with TOF [7,8].

The multifactorial etiology of TOF may contribute to disorders of transcription regulation networks [9,10]. TOF exhibiting genetic heterogeneity, including various forms such as single-nucleotide variants, copy number variations, transcriptional dysregulation, abnormal epigenetic modifications, and their complex interactions, all of which ultimately lead to gene network dysfunction [10,11,12,13]. Marcel and colleagues showed that deleterious private and rare mutations are in genes essential for the development of right ventricle and outflow tract in TOF patients, and that the interaction of these mutate genes form a similar abnormal gene regulatory network [12]. Similarly, Reuter et al. found that affected genes in TOF patients form an interactive network in which NOTCH1 plays a core role [13]. Weighted gene co-expression network analysis (WGCNA), a systems biology analysis method, has been widely used to reveal the gene modules and networks underlying abnormalities in various diseases [14,15,16]. Therefore, it is necessary to exploit systems biology analysis to synthesize the interactions between these genetic factors and their effects on TOF.

The maturity of the outflow tract (OFT) development that is closely related to TOF mainly depends on the neural crest cells (NCC) migrating towards the heart and the anterior second heart field cells (aSHFC) from the pharyngeal mesoderm, as well as involving various cell-cell interactions [17,18]. Profiling spatiotemporal gene expression of the NCC and aSHFC provides a new perspective for understanding the embryonic origin of TOF. Thus, the present study will further uncover the single-cell expression profile of Hub genes screened by WGCNA in the embryonic heart.

2 Materials and Methods

2.1 Data Availability and Preprocessing

All publicly available sequencing datasets related to TOF were retrieved from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/). Four datasets were ultimately included: two microarray expression profiles (accession numbers GSE26125 [19] and GSE35776 [20]), one transcriptomic dataset (GSE36761 [12]), and one methylation dataset (GSE62629 [21]). The preprocessed matrix data were downloaded from the GSE26125 and GSE35776, the original gene count expression matrix was downloaded from the GSE36761, and the original CpG score data were downloaded from the GSE62629. The analytical workflow is summarized in Fig. 1.

Figure 1: Flowchart of the data collection and analysis.

The GSE26215 data were normalized, and the GSE35776 expression matrix was subjected to log2 logarithmic transformation. The original expression matrix of GSE36761 was normalized by the vst function in the DESeq2 package (version 1.34.0; https://bioconductor.org/packages/release/bioc/html/DESeq2.html, accessed on 01 April 2025). Both GSE36761 and the methylation data were annotated using the GRCh37/hg19 genome assembly. To reduce technical noise, only genes exhibiting the top 60% intergroup variance between TOF and control samples in the GSE36761 expression matrix were retained for subsequent co-expression network construction.

2.2 Constructing Gene Network

A scale-free WGCNA network was constructed for the GSE36761 dataset using the WGCNA package (version 1.63) in R 4.0 (https://www.r-project.org, accessed on 01 April 2025), with the scale-free topology fit index exceeding 0.7. Detailed methodological protocols for WGCNA implementation can be found in our previous publication [22]. First, pairwise gene comparisons were performed to calculate absolute Pearson correlation coefficients. A soft threshold power β was determined as the minimal integer enabling scale-free network conformity. The adjacency matrix was then computed using the formula AM ab = cor a , b β , where AMab represents adjacency values between gene pairs a and b. We determined that β = 9 was the optimal weighting coefficient, i.e., the soft threshold. Subsequently, topological overlap matrices (TOM) were derived from the adjacency matrix to quantify network interconnectedness beyond direct correlations. Finally, hierarchical clustering with the average connectivity method was performed on the dissimilarity topological matrix to identify gene modules. The minimum number of genes in a module was set to 100, and similar genes were divided into the same module. Module eigengenes (ME), defined as the first principal components of module expression patterns, were calculated for dimensionality reduction. Highly correlated modules (r > 0.75) were merged using the dynamic tree cut algorithm.

2.3 Detecting TOF Modules

We assessed module-clinical trait associations by computing Pearson correlations between ME and clinical parameters (including age, gender, oxygen partial pressure, and disease status) to identify key modules of interest. The WGCNA package implements module significance (MS) as a metric quantifying module-trait relationships, calculated as the mean gene significance (GS) across all module genes. In our study, the module with maximal MS value for TOF status was designated as the TOF-associated module. Statistical significance was defined as p < 0.01 for MS values.

2.4 Enrichment Analysis

Pathway enrichment analysis was conducted using Metascape (v3.5.20220101; https://metascape.org, accessed on 01 April 2025) to annotate biological functions, incorporating Gene Ontology Biological Processes (GOBP), KEGG, Reactome, and WikiPathways databases. Gene Set Variation Analysis (GSVA) was performed with the GSVA package (v1.42.0) in R under default parameters. The GOBP gene sets were retrieved from the Molecular Signatures Database (MSigDB; http://www.gsea-msigdb.org/gsea/downloads.jsp, accessed on 01 April 2025). Differential pathway analysis between TOF and control groups was executed using the limma package (v3.50.1), with false discovery rate (FDR) correction applied for multiple testing. The simplifyEnrichment package (v1.4.0) subsequently clustered significant GOBP terms based on semantic similarity, generating condensed visualizations to facilitate rapid interpretation of predominant biological processes.

2.5 Identification of Hub Genes

To robustly identify Hub genes, we integrated protein-protein interaction (PPI) networks with WGCNA analysis through a dual-filter strategy. These filtered genes were subsequently used to construct a PPI network via STRING database (v11.5; https://string-db.org, accessed on 01 April 2025) with medium confidence threshold (combined score ≥ 0.4). Hub genes were defined as nodes with ≥5 interaction partners in the PPI network according to previous reports [23]. All network visualizations were generated using Cytoscape 3.9 (https://cytoscape.org, accessed on 01 April 2025), with node degree represented as topological centrality measure.

2.6 DNA Methylation and Single-Cell Transcriptomic Analysis

Differentially methylated CpGs analysis was performed using the R package methylKit (v1.20.0). Spatial scRNA-seq analysis was performed via the online database (EGAS00001003996) [24], which used to display the spatiotemporal expression information of the Hub genes in the embryonic heart. A murine single-cell dataset (GSE126128) encompassing E7.75-E9.5 cardiac development was analyzed to further validate Hub gene expression dynamics [25].

2.7 Data and Code Availability

All datasets analyzed in this study were obtained from the aforementioned repositories. The WGCNA network construction workflow was implemented using code adapted from the official WGCNA tutorial (Version 1.72-1; https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/, accessed on 01 April 2025), with full reproducibility ensured through open access to all analytical scripts.

3 Results

3.1 WGCNA Network Construction

After hierarchical clustering analysis, TOF-18 was identified as an outlier sample subsequently excluded to mitigate confounding effects. As shown in Fig. 2, the reclustered samples demonstrated distinct segregation between the case and control groups. The adjacency matrix transformation was performed on the 10,503 genes of the 28 samples, and then the transformation into a topological matrix was continued (Fig. 3A–C), confirming the resultant network exhibited scale-free topology characteristics. In addition, module detection was implemented using dynamic tree cutting with a minimum module size of 100 genes, followed by merging of highly correlated modules (similarity threshold > 0.75). This process ultimately yielded eight distinct co-expression modules (Fig. 3D).

Figure 2: Sample clustering tree and clinical information heatmap. The branches in the clustering tree correspond to the samples in the squares in the heatmap below.

Figure 3: Soft threshold and clustering tree diagram of the gene module. (A) Relationship analysis between different soft thresholds (weighting coefficients) and network topology. (B) Relationships between different soft thresholds and the average connectivity in the network. (C) Correlation scatter plot of log10(k) and log10(p(k)). (D) Clustering tree diagram of all genes.

3.2 Module-Trait Association Analysis

We systematically evaluated the associations of eight co-expression modules with clinical characteristics (age, sex, oxygen saturation, and disease status) to identify clinically relevant modules. Notably, the ME of the darkorange and blue modules demonstrated highly significant correlations with TOF, with correlation coefficients of 0.84 and 0.78, respectively (Fig. 4A). In addition, the darkorange module also presented the highest MS value (Fig. 4B). Although the blue module ranked second in MS and showed statistically significant TOF correlation, its intramodular connectivity metrics within the WGCNA framework were suboptimal. Consequently, only the darkorange module was retained as the clinically meaningful, TOF-associated module for downstream analyses.

Figure 4: Identification of clinically significant gene modules. (A) Module-trait association heatmap demonstrating correlations between clinical features (x-axis) and module eigengenes (ME) values (y-axis). (B) Module significance histogram indicating the average gene significance (GS) values across eight modules.

3.3 Identification of Hub Genes in the Darkorange Module

Hub gene identification within the darkorange module was performed by analyzing intramodular connectivity metrics. For the 481 genes assigned to this key module, a robust correlation was observed between GS and module membership (MM), reinforcing the biological relevance of this module to TOF pathogenesis (Fig. 5A). Functional enrichment analysis via Metascape revealed predominant associations with organogenesis (e.g., cardiac, renal, and cranial development; epithelial morphogenesis; cell junction assembly) and WNT signaling-mediated cardiac development, alongside system-level processes (cardiac contraction and circulatory regulation) (Table 1). PPI network analysis using the STRING database identified eight Hub genes (GJA1, PTK7, PRICKLE1, SFRP1, RUNX2, FBN2, PIK3R1, and KCNJ2), each exhibiting ≥5 direct interactions within the network (Fig. 5B).

Figure 5: Identification of Hub genes. (A) Scatter plot demonstrating the correlation between module membership (MM) and gene significance (GS) values. (B) Protein-protein interaction (PPI) network of 185 WGCNA-identified Hub genes analyzed via STRING database.

3.4 Dysregulated Cell Polarity-Related Signaling Pathway in TOF

After the identification of the Hub genes, we observed that five Hub genes were significantly enriched in pathways governing epithelial morphogenesis, in which SFRP1, PRICKLE1 and PTK7 were mainly involved in the establishment of planar cell polarity in epithelial morphogenesis (Table 1). To further validate these results, we used two TOF transcriptomic datasets (GSE26125 and GSE35776) for GSVA. The results suggest that epithelial polarization-related pathways (embryonic epithelial morphogenesis pathway and anterior-posterior axis specification in embryos pathway) and their upstream noncanonical WNT signaling are also inactivated, whereas the canonical WNT signaling are activated (Fig. 6A). Therefore, the next analysis focused on five genes involved in the morphological and polarity development of embryonic epithelial cells: GJA1, PTK7, PRICKLE1, SFRP1 and RUNX2. Clustering of 1357 dysregulated biological processes identified through GSVA revealed predominant involvement in metabolic regulation, organogenesis, signal transduction, macromolecular assembly in organelles, and cell cycle/proliferation (Fig. 6B). These findings align with established pathobiological mechanisms underlying congenital cardiac malformations, reinforcing the hypothesis that disrupted epithelial polarity and WNT signaling dynamics contribute to TOF pathogenesis.

Figure 6: Differential pathway activity and functional clustering between TOF and healthy control. (A) Differential pathway activity scores between TOF and controls quantified by Gene Set Variation Analysis (GSVA) across GSE26125 and GSE35776 datasets. (B) Functional similarity clustering of differentially active pathways from GSVA.

3.5 Methylation Modifications on GJA1, PRICKLE1, PTK7, and RUNX2 in TOF

To investigate the potential causes of epithelial morphogenesis pathway disruption, we analyzed the mutational landscapes of GJA1, PTK7, PRICKLE1, SFRP1 and RUNX2 in two whole-exome sequencing cohorts (the US CHD cohort [26] and the UK TOF cohort [27]). In the US CHD cohort, two novel de novo variants were identified in conotruncal defect patients: a missense mutation (GJA1:c.121768138C>G) and a synonymous variant (PTK7:c.43106663C>T). Additionally, a loss-of-function mutation in PTK7 and SFRP1 was detected in control subjects. The UK cohort revealed two rare RUNX2 mutations (Table 2). Given the rarity of coding sequence alterations, we hypothesized alternative regulatory mechanisms and interrogated methylation profiles using the GSE62629 dataset. Strikingly, promoter hypermethylation was consistently observed in TOF cases for GJA1, PRICKLE1, RUNX2, and PTK7, suggesting transcriptional repression of these polarity regulators (Table 3).

3.6 Expression of GJA1 and SFRP1 in Embryonic OFT-Derived Cells

Systematic analysis of the candidate genes in TOF bulk RNA-seq data (GSE36761, 22 TOF cases vs 8 controls) revealed: All candidate genes (GJA1, PTK7, PRICKLE1, SFRP1 and RUNX2) showed significant upregulation in TOF myocardial tissues, as visualized in Supplementary Fig. S1.

We used two single-cell datasets to analyze the expression of these genes in different locations and at different times in the embryo. The results revealed the selective expression of GJA1 and SFRP1 in human embryonic hearts at 4.5–9 gestational weeks (Carnegie CS13-CS23) (Table 4) [24]. GJA1 was expressed in trabecular ventricular myocardium, whereas SFPR1 was expressed in OFT/large vessels, mediastinal mesenchyme and vessels (Table 4). Similar expression patterns were also found in mouse hearts at the embryonic age of E7.75–9.25 days (Carnegie CS9-CS13) [25]. GJA1 demonstrated ubiquitous expression across critical cardiac lineages, including multipotent progenitors, NCC, endocardial cells, lateral plate mesoderm (LPM), and paraxial mesoderm (PM) cells (Fig. 7). SFRP1 exhibited stage-specific enrichment, predominantly in multipotent progenitors, mesodermal structural cells, and NCC, LPM and epicardial cells at E8.25 and E9.25 (Fig. 6). PTK7 was only slightly expressed in the LPM and PM, RUNX2 was virtually absent in the heart, and PRICKLE1 was expressed in trace amounts in the heart (Fig. 7).

Given the essential roles of NCC and aSHFC in the development of the OFT, we observed the expression patterns of five Hub genes. Our analysis revealed that GJA1 was significantly expressed in the aSHFC-derived multipotent progenitors, posterior second heart field (pSHF) cells, and branchiomeric muscle progenitors, with progressive downregulation culminating in near-complete silencing at E9.5 (Fig. 7). In contrast, SFRP1 maintained stable expression in aSHF and pSHF, and both SFRP1 and GJA1 began to be expressed in large quantities in NCC at E8.25 (Fig. 7). PTK7 and PRICKLE1 were weakly expressed in these three types of cells, whereas RUNX2 was not expressed (Fig. 8). Thus, the GJA1, SFRP1, PRICKLE1 and PTK7 genes were ultimately identified as candidate pathogenic genes of TOF.

Figure 7: Expression profiles of Hub genes in the hearts of mice at E7.57, E8.25, and E9.25.

Figure 8: Expression profiles of Hub genes in mouse hearts at embryonic ages E7.57, E8.25, and E9.25.

4 Discussion

In this study, we constructed WGCNA network to identify disease-associated modules, subsequently prioritizing eight Hub genes through PPI network analysis. Functional enrichment revealed five genes (GJA1, PTK7, PRICKLE1, SFRP1 and RUNX2) converging on epithelial morphogenesis pathways, particularly planar cell polarity (PCP) signaling. Epigenetic profiling demonstrated promoter hypermethylation in GJA1, PRICKLE1, PTK7 and RUNX2, aligning with observed PCP pathway inactivation. In addition, developmental validation across NCC and aSHFC lineages confirmed the GJA1 and SFRP1 expression, while RUNX2 showed no detectable expression in human or murine cardiac progenitors. Thus, the GJA1, SFRP1, PRICKLE1 and PTK7 genes were ultimately identified as candidate pathogenic genes of TOF. Our findings implicate epithelial polarity pathways as central mechanistic drivers of TOF, highlighting their therapeutic potential for congenital heart disease.

Cardiac development initiates from a linear epithelial tube with arteriovenous polarity, requiring precisely orchestrated epithelial morphogenesis involving cellular polarity establishment. Cell polarity is a basic attribute of various cells and manifests as asymmetry in morphology, macromolecule distribution and function, including apicobasal polarity (ABP) and PCP, which play important roles in embryonic organogenesis and development [28,29,30,31]. Notably, aSHFC exhibit hallmarks of an atypical ABP epithelium, with a single cilium at the apex and basally enriched dynamic actin-rich filopodia [32,33]. Disruption of these polarized features triggers aSHFC proliferative arrest and ectopic differentiation, ultimately impairing heart tube elongation and OFT morphogenesis [32,34,35,36,37]. In addition, epithelial cell polarity dysregulation has been shown to induce embryonic OFT truncation, perturbing ventriculoarterial alignment and culminating in the cardinal pathological features of TOF [38,39]. Therefore, the dysregulation of epithelial cell polarity and related pathways is a potential mechanism of the embryopathogenesis of TOF.

GJA1 encodes connexin 43 (Cx43), which is a critical component of myocardial gap junctions that mediates intercellular communication [40]. Online Mendelian inheritance in humans catalogs pathogenic information on GJA1 (MIM* 121014), which associated with ventricular septal defects and hypoplastic left heart syndrome. A cohort study of 152 TOF patients identified GJA1 point mutations in only 5.3% of cases (8/152). however, murine models showed neither congenital heart defects nor reduced viability upon monoallelic or biallelic Gja1 knockout [41]. Paradoxically, Rhee et al. reported that Gja1 knockout mice exhibit conus malformation and coronary artery abnormalities, with mechanistic analyses revealing Cx43-microtubule interactions as essential regulators of cellular polarity [42]. Another study has demonstrated that cardiomyocytes in children with TOF exhibit reduced Cx43 expression and an abnormally diffuse distribution across the cell surface, resembling patterns seen in immature tissues [43]. Thus, these observations suggest potential genetic mechanisms underlying GJA1 dysregulation that extend beyond conventional genetic sequence variations, possibly involving DNA methylation patterns identified in our investigation or alternative epigenetic regulatory pathways.

SFRP1, a soluble member of the SFRP family, functions as a canonical WNT signaling antagonist by competitively binding to WNT ligands through its conserved frizzled-like cysteine-rich domain [44,45]. During murine cardiac development (between E8.5 and E12.5), SFRP1 and WNT8 are localized specifically in cardiomyocytes, but neither detected in the pericardium nor the endocardium [46]. This localization pattern aligns with our observations in human and mouse embryonic hearts within the cardiogenic zone and OFT. Co-immunoprecipitation assays confirmed direct binding between SFRP1 and WNT8, establishing their functional interplay during cardiogenesis [46]. Gibb et al. reported that differentiating cardiomyocytes robustly upregulate SFRP1 expression to antagonize canonical WNT signaling, thereby facilitating OFT morphogenesis [47]. These findings suggest a critical role for SFRP1 in OFT development and its potential involvement in TOF pathogenesis.

The PRICKLE1 protein, a core component of the PCP pathway, is essential for establishing both PCP and ABP in embryonic cells [48]. PRICKLE1 knockout mice exhibit failures in distal visceral endoderm migration and mesoderm formation, leading to early embryonic lethality [48]. Additional studies demonstrate that embryonic loss of PRICKLE1 disrupts cellular polarity, causing disorganized myocardial fiber alignment and OFT malformations [49,50]. Chiapparo et al. further associated PRICKLE1 deficiency with OFT defects, identifying it as a downstream target of MESP1 that regulates migration velocity, polarity, and directional movement of cardiovascular progenitor cells [51].

PTK7, a member of the protein tyrosine kinase (PTK) family, is a transmembrane glycoprotein critical for embryonic development and tissue homeostasis [52,53]. Structurally, it comprises seven extracellular immunoglobulin-like domains, a transmembrane region, and an intracellular pseudokinase domain lacking catalytic activity, enabling its role as a co-receptor and signaling scaffold [52]. As a core PCP regulator, PTK7 coordinates tissue morphogenesis via non-canonical WNT signaling [54]. Studies in mice reveal that Ptk7 mutations can induce severe neural tube defects alongside multisystem anomalies, including cardiac malformations such as double outlet right ventricle and ventricular septal defects [55,56]. However, the role of PTK7 in OFT development and TOF pathogenesis requires further investigation.

This study has several limitations. First, as a retrospective investigation, its conclusions are inherently constrained by the sample size limitations of the database analyzed; future prospective studies with larger, multi-center cohorts are warranted to validate these findings. Second, while DNA methylation is a well-studied epigenetic modification, other mechanisms such as histone modifications, non-coding RNA regulation, and chromatin accessibility also play critical roles in gene regulation. Thus, it is necessary to explore these additional layers of epigenetic regulation in future studies. Finally, while partial molecular mechanisms of GJA1 and PRICKLE1 have been documented in prior research, comprehensive validation through integrated in vitro and in vivo experimental models remains essential to elucidate the detailed pathogenic mechanisms and functional interplay among all four candidate genes identified in this work.

5 Conclusion

In summary, our study systematically dissects dysregulated signaling networks in in TOF. We identified GJA1, SFRP1, PRICKLE1, and PTK7 as critical modulators of epithelial cell polarity in TOF. These coordinated epigenetic, transcriptional, and functional perturbations collectively define a pathogenic cascade that translates cellular polarity defects into macroscopic cardiac structural malformations. Our findings implicate epithelial polarity pathways as central mechanistic drivers of TOF, highlighting their therapeutic potential for congenital heart disease.