FASTQ files for all probands and their parents were assessed for quality control using FASTQC and MultiQC (Figure 1A) (39, 40). 90% of the FASTQ files for the selected TOF trios and 91% of SFARI control trios passed quality control (i.e., the proband and both parents passed) and were used for our analyses of de novo variants (PCGC TOF n = 362; SFARI control n = 317). Files were aligned with BWA-MEM to GRCh38 and then processed according to GATK Best Practices workflow (41, 42). Variants were called using GATK HaplotypeCaller. Variants were filtered for only high confidence de novo variants using GATK CalculateGenotypePosteriors, wherein both parents had a confident call for being homozygous for the reference allele. Functional annotation was next performed with snpEff and filtered with snpSift (43, 44). Additional functional annotation was conducted with ANNOVAR, including for MAF annotation derived from the Exome Aggregation Consortium (ExAC) and Combined Annotation-Dependent Depletion (CADD) scoring (45–47). CADD scores can score SNPs and indels based on their likelihood to be causal both on factors such as species conservation and impact on protein functionality. This analysis pipeline was automated utilizing Snakemake (48, 49). The complete bioinformatic pipeline for this analysis can be accessed via github¹.

Variants annotated as synonymous, upstream, downstream, 5′ or 3′ UTR, and intronic were filtered out (Figures 1A,B). Variants in genes known to be false positive calls in exome sequencing research were also filtered out (50). To focus on rare alleles exclusively we selected a MAF threshold of 1 × 10⁻⁴ according to ExAC via annotation by ANNOVAR. Only variants annotated by snpEff to be of moderate or high impact and a defined CADD phred score of 20 (representing approximately the top 1% of single nucleotide variants that had the highest prediction of functional effect) or greater were included in the enrichment analysis, with a secondary analysis conducted on all variants with a moderate or high impact annotation (snpEff) but no inclusion cutoff based on CADD phred scores.

A combination of Gene Ontology (GO) tools was used to determine enrichment of biological processes and pathways: the Protein ANalysis THrough Evolutionary Relationships (PANTHER) classification system², the Database for Annotation, Visualization, and Integrated Discovery (DAVID) software³, Metascape⁴) and Qiagen's Ingenuity Pathway Analysis (IPA) software (51–55). The input to each of these tools was the variant list generated from our pipeline, one list of all de novo, non-synonymous variants and one list of only those with a defined CADD score of 20 or greater (Figure 1B). Genes identified in the variant analysis were evaluated for expression levels in the fetal heart via the Brotman Baty Institute's Atlas of Gene Expression During Development⁵ and those related to heart development were identified using IPA (56). All candidate variants were validated by visualization using Integrative Genomics Viewer (IGV) and for 65 probands that additionally had genome sequencing data available, variant calls from the exome sequence data were compared against variant calls from genome sequence data previously conducted by PCGC. None of the variants identified in the PCGC cohort was identified in the SFARI cohort, further supporting that effects seen are highly unlikely to be due to random enrichment of nonpathogenic alleles. Statistical analyses were conducted either by the GO analysis tools (PANTHER, DAVID, Metascape, IPA) or using the stats package in R (Version 4.0.2) (57).

The data and study materials are under restricted access, as the data use agreements with PCGC and SFARI prohibit distribution of patient-level data, but upon application may be available to other researchers for the purposes of reproducibility or replicating the procedure. Data can be requested from both PCGC⁶ and SFARI⁷ by researchers. Exome sequencing data from PCGC was accessed via dbGaP authorized access of study phs001194.v3.p2 and exome sequencing data from SFARI was collected from the Simons Foundation Powering Autism Research for Knowledge (SPARK) project.

Of the 362 patients in the final TOF cohort more individuals were assigned male at birth (n = 203) than female (n = 159; 56.08% vs. 43.92%). As these are just the cases of TOF that had exome sequencing data for the proband and both parents and the associated FASTQ files for all three passed initial quality metrics, this does not mean there is such a discrepancy amongst the total non-syndromic TOF population in the PCGC database. When looking at all non-syndromic conotruncal defects from PCGC however, a similar proportion was described with 57.6% described as male and 42.4% female (58). A similar breakdown was achieved in the SFARI cohort of 317 unaffected siblings of children with autism (53% male vs. 47% female). The racial and ethnic composition of the TOF cohort (Supplementary Figure 1A) is skewed more heavily white (83.1% vs. 72%), with the heaviest differences coming from far fewer Black individuals (5.8% vs. 13%) compared to the U.S. population at the time of sample collection and no individuals in our TOF cohort identified as American Indian, Alaska Native, Native Hawaiian or Pacific Islander. Unlike participants in PCGC's database, those in the SPARK project for SFARI did not have race or ethnicity data collected. Instead HapMap haplotypes were utilized to categorize ancestral background. The SFARI control cohort ancestral superclass composition was 3.42% African, 9.63% Admixed American, 1.24% East Asian, 79.81% European, 1.86% South Asian and 4.04% unknown ancestral superclass. While one cannot presume that the 79.81% EUR population is entirely white, the probability is high that the control cohort like the TOF cohort is predominately white.

The de novo mutation rate in the TOF cohort was 0.33 per proband for variants with a CADD score of 20 or greater (120 variants, 117 genes, 105 probands) and 0.56 per proband for all variants annotated as moderate or high impact by snpEff (183 variants, 179 genes, 147 probands; Supplementary Table 1). All variants identified after filtering are listed in Supplementary Table 1, with those with the 15 highest CADD scores listed in Table 1. The de novo mutation rate in our control cohort was 0.31 (97 variants, 97 genes, 82 probands) for those with a CADD score of 20 or greater and 0.51 overall (161 variants, 161 genes, 123 probands), neither of which was statistically different from the TOF cohort (Chi-square p = 0.88 and p = 0.47 respectively).

Examining variants with a CADD score of 20 or greater in the PCGC TOF probands, 11.67% were LoF mutations and 88.33% D-Miss/Indel, compared to all variants where 14.75% were LoF and 85.25% D-Miss/Indel (Supplementary Table 2). By comparison, there were significantly fewer LoF variants with a CADD score of 20 or greater in the SFARI control cohort (Supplementary Figure 1C 4.12%, 4/97, OR 3.07, 95% CI 0.98–9.66, p = 0.05). There was no statistical difference in LoF variants between the two cohorts when examining all variants with moderate/high impact from snpEff alone (13.04%, 21/161, OR 1.15, 95% CI 0.62–2.13, p = 0.65). An evaluation of gene expression of candidate genes using the Brotman Baty Institute's Atlas demonstrated that more than half of all genes affected both in the CADD 20 analysis (50.43%, 59/117) and analysis of all variants (53.63%, 96/179) were found to be expressed in more than 1% of fetal cardiac tissue (56). This is significantly more than the number of genes affected in the control cohort, where only slightly over a third of genes were found to be expressed in the fetal heart in both analyses (CADD 20 analysis: 37.11%, 36/97, OR 1.72, 95% CI 0.99–2.98, p = 0.05; overall analysis: 36.02%, 58/161, OR 2.05, 95% CI 1.33–3.18, p = 0.001).

GO analysis found CHD (7-fold enrichment, p = 5.01 × 10⁻⁷, Benjamini-Hochberg (B–H) p = 0.003) to be enriched in our analysis, including TOF (8.5-fold enrichment, p = 1.26 × 10⁻⁶, B–H p = 0.004) and CHARGE syndrome (a syndromic form of TOF; 21-fold enrichment, p = 3.98 × 10⁻⁵, B–H p = 0.04; Figure 2, Table 2). Only CHD was significant in the analysis of all variants (5-fold enrichment, p = 6.31 × 10⁻⁶, B–H p = 0.02). According to IPA Disease and Function analysis there was also a significant overlap with genes involved in morphogenesis of the cardiovascular system (p = 1.94 × 10⁻⁶, B–H p = 5.69 × 10⁻⁵) and congenital anomaly of the cardiovascular system (p = 1.1 × 10⁻⁶, B–H p = 3.03 × 10⁻⁵). Overlap with these genes was less significant in the analysis of all variants (morphogenesis p = 8.57 × 10⁻⁵, B–H p = 1.67 × 10⁻³; anomaly p = 3.39 × 10⁻⁴, B–H = 5.85 × 10⁻³). The control cohort did not have overlap with any of these terms related to CHD or cardiac development.

A total of four genes were mutated in two probands in the TOF cohort: CHD7, CLUH, UNC13C, and WASHC5 (Table 3). CHD7 is a helicase critical in chromatin remodeling, CLUH is a mRNA binding protein involved in mitochondrial regulation, UNC13C is a membrane fusion protein involved in vesicle maturation and synaptic transmission and WASHC5 is a component of the WASH complex that functions in endosomes (59–62). Variants in these genes were annotated as having a high impact by snpEff or an average CADD score greater than 20. No variant was identified in more than two probands. Of the genes mutated in two probands, WASHC5, CLUH and CHD7 were identified to be expressed in the heart, of which WASCH5 and CHD7 have previously been identified to be related to CHD pathogenesis, including syndromic cases of TOF (8, 11, 14, 31, 63–65). Our findings further the case that these genes are also important in non-syndromic TOF. CLUH and UNC13C represent novel candidate genes that have not been previously associated with TOF pathogenesis, but are of high interest due to their observed perturbation in more than one unrelated probands.

The rate of likely damaging de novo variants (0.56) in our study is comparable to other studies of rare, de novo variants in cohorts of children with CHD (31, 66). Our analysis also identified variants in genes, pathways and biological processes previously identified in other studies of TOF. CHD7, JAG1, and WASHC5 have been implicated in TOF, and ATP2A2, CHD4, DNAH6, DTNA, EDNRA, FAT4, MED13l, RASA1, and SLC5A2 in CHD generally according to the literature and our IPA analysis (Supplementary Table 3) (8, 10, 11, 14, 35, 63, 64, 66, 67). A number of these genes have been established to be critical to the progenitor populations necessary for heart development: first and second heart field derived cells and cardiac neural crest cells. EDNRA is expressed for instance in the first heart field and EDNRA null mice have hypoplastic ventricles (68). MED13l is expressed in migratory cardiac neural crest cells, is heavily associated with transposition of the great arteries and coarctation and in cardiac energy metabolism and apoptosis (69). FAT4 is expressed in cardiac mesoderm that produces both the first and second heart fields and has been demonstrated to regulate cardiac progenitor proliferation and planar cell polarity (70). Four genes had de novo variants that were found in multiple probands in our cohort, of which two (CHD7 and WASHC5) are in genes previously implicated in TOF (Table 3). This finding speaks both to the centrality of these genes in TOF biology and the ability of our analysis to detect variants with high likelihood for pathogenicity.

Many exome sequencing studies of TOF have identified an enrichment for variants in genes related to autism, chromatin remodeling, cilia and ciliopathy, VEGF signaling, Notch signaling, and the Wnt/β-catenin pathway (71–77). Chromatin remodeler activity (19.61-fold enrichment, p = 7.54 × 10⁻⁵, B–H p = 0.025) was indeed found to be enriched in our cohort (Table 2), as was a significant relationship found for genes related to autism spectrum disorder and intellectual disability (p = 2.18 × 10⁻⁶, B–H p = 5.74 × 10⁻⁵). It is worthy to note that autism spectrum disorder and intellectual disability (p = 0.003, B–H p = 0.03) was also found to be enriched in our control cohort though to a lesser degree. This may be in relation to the large number of genes encompassed in this disease (SFARI's Gene database lists 942 genes as of 2022) and a reflection of the disagreement between datasets and research groups on what true risk genes are for autism spectrum disorder and intellectual disability (78–80). Ciliary body morphogenesis (GO:0061073; 42.89-fold enrichment, p = 0.03, B–H p = 1.00), the Notch signaling pathway (P00045; 7.63-fold enrichment, p = 0.03, B–H p = 0.81) and the Wnt signaling pathway (P00057; 3.82-fold enrichment, p = 0.003, B–H p = 0.43) were found to be enriched in the variants called in our TOF cohort; however, these pathways were not significant after correction for multiple comparison testing.

Molecular function GO analysis demonstrated that ATP binding (2.75-fold enrichment, p = 5.61 × 10⁻⁶, B–H p = 0.005) and DNA helicase activity (11.59-fold enrichment, p = 9.62 × 10⁻⁵, B–H p = 0.03) were found to be enriched in our identified variants with a CADD score of 20 or greater (Figure 2, Table 2). Both terms were also significantly enriched in the analysis of all variants, in addition to meiotic cell cycle (GO:0051321; 6.83-fold enrichment, p = 2.40 × 10⁻⁶, B–H p = 0.05). Biological function GO analysis elicited that aorta development (21.9-fold enrichment, p = 3.39 × 10⁻⁶, B–H p = 0.04), cardiac chamber development (10.9-fold enrichment, p = 4.00 × 10⁻⁶, B–H p = 0.04), and general cell development (2.78-fold enrichment, p = 6.19 × 10⁻⁷, B–H p = 0.01) were significantly enriched. By comparison, no functional GO terms were significantly enriched in the analysis of all variants. Diseases identified to be enriched related to our variants are myopathy (5.7-fold enrichment, p = 1.58 × 10⁻⁷, B–H p = 0.001), congenital ear anomaly (13-fold enrichment, p = 1.00 × 10⁻⁶, B–H p = 0.004) and developmental disabilities (7.3-fold enrichment, p = 1.26 × 10⁻⁶, B–H p = 0.004).

IPA analysis detected a significant overlap with genes involved in cell development, neural development, muscle development and function, lipid metabolism, hearing impairment (p = 5.93 × 10⁻⁶, B–H p = 1.42 × 10⁻⁴), perinatal death (p = 6.53 × 10⁻⁵, B–H p = 0.001), development of the body axis (p = 1.16 × 10⁻⁴, B–H p = 0.002), morphogenesis of semilunar valve (p = 5.3 × 10⁻⁴, B–H p = 0.008), and development of the heart septum (p = 0.003, B–H p = 0.03). Processes involved in cell development include assembly of the spindle apparatus (p = 7.13 × 10⁻⁵, B–H p = 0.001), adhesion of cell-associated matrix (p = 1.87 × 10⁻⁴, B–H p = 0.003), cell cycle progression (p = 9.81 × 10⁻⁴, B–H p = 0.013), and organization of cytoskeleton (p = 0.001, B–H p = 0.014). Included in neurocognitive development are familial central nervous system disease (p = 6.14 × 10⁻⁶, B–H p = 1.46 × 10⁻⁴), proliferation of neural cells (p = 0.001, B–H p = 0.015), innervation of neurons (p = 0.001, B–H p = 0.016), growth of axons (p = 0.0015, B–H p = 0.018), and guidance of axons (p = 0.003, B–H p = 0.034). The analysis of all variants revealed an additional enrichment for differentiation of neural cells (p = 6.95 × 10⁻⁴, B–H p = 0.05). Related to muscle development and function are transport of calcium (p = 3.76 × 10⁻⁵, B–H p = 7.53 × 10⁻⁴), organization of muscle cells (p = 4.09 × 10⁻⁴, B–H p = 0.006), and muscle contraction (p = 0.001, B–H p = 0.017). We have grouped calcium transport into the category of muscle development and function as we believe that it is distinctly tied to dysfunction of the sarcoplasmic reticulum. ATP2A2 and RYR1 are shared between calcium transport, muscle contraction and organization of muscle, with ANXA6 and EDNRA additionally shared between calcium transport and muscle contraction. While not significant after multiple comparison testing, both calcium regulation in cardiac cells (WP536; 8.49-fold enrichment, p = 3.24 × 10^–4, B–H p = 0.25) and sarcoplasmic reticulum calcium ion transport (GO:0070296; 22.88-fold enrichment, p = 0.004, B–H p = 0.65) specifically were enriched. Related to lipid metabolism, quantity of long-chain acyl-coenzyme A (p = 8.57 × 10⁻⁴, B–H p = 0.01) and metabolism of ceramide (p = 9.52 × 10⁻⁴, B–H p = 0.01) were found to be significant, with the additional terms of accumulation of phosphatidylinositol 3,4-diphosphate (p = 0.001, B–H p = 0.02), accumulation of phospholipid (p = 0.002, B–H p = 0.02), and accumulation of phosphoinositide (p = 0.004, B–H p = 0.05) significantly enriched in the analysis of all variants. A majority of the enrichments from the CADD 20 analysis were found in the IPA analysis of all variants (Supplementary Table 3). Of note, terms associated with calcium handling and muscle contraction were only identified in the CADD 20 analysis (concentration of Ca⁺², liberation of Ca⁺², oscillation of Ca⁺², transmembrane transport of Ca⁺², formation of muscle, hereditary myopathy, and muscle contraction; p-values and B–H values in Supplementary Table 3), as was development of heart septum, demonstrating the critical nature of the CADD 20 analysis.

Of final note, there was also significant enrichment of target genes of four different transcription factors, all of which have previously been tied to cardiac development or injury response. The most significant is SOX9 (7.6-fold enrichment, p = 3.98 × 10⁻⁵, B–H p = 0.04), which coincidentally is the only transcription factor that has been previously reported in relation to TOF and is critical for ventricular septum formation (81, 82). NKX6.2 (7.3-fold enrichment, p = 5.01 × 10⁻⁵, B–H p = 0.05) has been previously tied to CHD. While no specific role of NKX6.2 in cardiac development has been identified, the Nkx family (particularly NKX2.5) plays a critical role in specification of the second heart field pharyngeal mesoderm (83, 84). HFH1 (7.2-fold enrichment, p = 6.31 × 10⁻⁵, B–H p = 0.05), also known as FOXQ1, is involved in cell proliferation, differentiation, and myocardial fibrosis, though it has mainly been studied in cancer (85, 86). HFH4 (6.9-fold enrichment, p = 7.94 × 10⁻⁵, B–H p = 0.05) has been implicated in ciliogenesis. This process by itself is not enriched in our TOF cohort but has been reported previously in other genomic analyses of TOF and left-right body axis formation, which we demonstrated is enriched in our cohort (34, 87). In contrast, none of the four transcription factors enriched in the CADD 20 variant list is significant after multiple comparison correction in the overall list. This again demonstrates the importance of the CADD 20 analysis in identifying the variants and their associated pathways that would not be elucidated from the analysis of all de novo non-synonymous variants in the cohort. IRF is the only transcription factor of significance (6.2-fold enrichment, p = 1.58 × 10⁻⁵, B–H p = 0.03). IRF is enriched in the CADD 20 analysis (5.3-fold enrichment, p = 0.003, B–H p = 0.46), though not after correction for multiple comparison testing. IRF is predominately connected to cardiac biology through its role in cardiac fibrosis, ventricular remodeling, and heart failure (88, 89, 90). The fact it is significantly enriched in the larger analysis, that it has a functional similarity to HFH1/4, and the stringency of our various filtering metrics warrants further consideration of its role in TOF pathogenesis and prognosis.

Enrichment analyses of the control cohort resulted in no identified enriched pathways outside of analysis with IPA. The only overlap in the IPA analysis was for enrichment in growth of axons (p = 5.47 × 10⁻⁴, B–H p = 0.009), proliferation of neural cells (p = 0.002, B–H p = 0.02), organization of cytoskeleton (p = 0.005, B–H p = 0.04), and autism spectrum disorder or intellectual disability (p = 0.003, B–H p = 0.03). While growth of axons was more significantly enriched than in the TOF cohort, all other terms were enriched to a lesser degree than the TOF cohort. Whether or not these should not be considered as potentially implicated in TOF pathogenesis requires further study.