Abstract

   Tetralogy of Fallot (TOF), the most common severe cyanotic congenital
   heart disease, has unclear genetic causes. Through next-generation
   sequencing in 131 patients with nonsyndromic TOF, we identified an
   increased burden of rare deleterious variants in ciliary genes and
   cilium pathway and observed a multigenic inheritance pattern, with an
   odds ratio (OR) of 1.672 [95% confidence interval (CI), 1.120 to 2.547;
   P = 0.0104] for more than two deleterious variants and a cumulative OR
   of 3.158 (95% CI, 1.381 to 6.371; P = 0.0038) for six variants.
   Functional validation in single- and double-heterozygous mouse models
   carrying these variants recapitulated TOF-like phenotypes and impaired
   normal cilia structure and function, particularly disrupting Hedgehog
   signaling in cardiomyocytes, and down-regulating key transcription
   factors Gata4 and Nkx2-5. Together, our study provides compelling
   evidence linking ciliary gene variants to a heightened risk of TOF in
   Han Chinese, offering valuable genetic insights into the etiology and
   pathogenesis of nonsyndromic TOF and supporting a multigenic
   inheritance model for the disease.
     __________________________________________________________________

   Ciliary gene variants contribute to TOF risk, revealing a multigenic
   inheritance pattern and insights into TOF pathogenesis.

INTRODUCTION

   Congenital heart disease (CHD) is the most common type of congenital
   malformation, accounting for 28% of all major congenital anomalies
   ([62]1) and affecting 8 per 1000 live births ([63]2). Tetralogy of
   Fallot (TOF) is the most frequent cyanotic CHD, with a global
   prevalence of about 1 in 3000 live births ([64]3). TOF is considered a
   malformation of the cardiac outflow tract (OFT) that includes
   structural abnormalities, i.e., overriding aorta (OA), ventricular
   septal defects (VSD), and pulmonary artery (PA) stenosis, and
   compensated right ventricular hypertrophy ([65]4). With recent advances
   in the diagnosis, surgical treatment, and postoperative care, almost
   all of those born with TOF can now be expected to survive to adulthood
   after successful surgical correction at a young age ([66]4).
   Unfortunately, as long-term sequelae, patients with TOF have high
   morbidity of arrhythmia and ventricular dysfunction after 40 years of
   age ([67]5, [68]6), which seriously affect quality of life and bring a
   heavy burden to the society.

   TOF is a complicated disease influenced by both genetic and
   environmental factors ([69]7). Approximately 20% of patients with TOF
   present with additional extracardiac anomalies and are classified as
   syndromic TOF ([70]8). Among them, 22q11.2 deletion syndrome (DiGeorge
   syndrome) and trisomy 21 (Down syndrome) account for approximately 15
   and 7% of cases, respectively ([71]9). Multiple pathogenic loci, such
   as chromosomal copy number variations (CNVs) at 22q11.2 and mutations
   in genes like GATA4, have been identified in syndromic TOF and provided
   valuable insights into its underlying mechanisms ([72]9, [73]10).
   Nevertheless, nearly 80% of cases are sporadic and nonsyndromic TOF,
   with no clearly identified etiology or major mutated genes ([74]11,
   [75]12). A recent large-scale study involving nine research centers
   across Europe and Australia identified NOTCH1 as the most frequently
   mutated genes in nonsyndromic TOF, followed by FLT4 ([76]13). To date,
   despite the well-established function of several genes in the
   pathogenesis of TOF, only a small fraction of nonsyndromic TOF cases
   can be attributed to a detectable abnormality in any one of the known
   TOF genes identified: CNVs in the TBX1 gene located in the 22q11.2
   region (approximately 6%) ([77]14) and NKX2.5 mutations (approximately
   4%) ([78]15). The underlying causes remain unknown for most patients
   with nonsyndromic TOF. In addition, the lack of experimental profiling
   to determine the biological mechanisms involved in these genetic
   variants hampers our in-depth knowledge of the disease.

   To systematically investigate the genetic etiology of nonsyndromic TOF,
   we performed next-generation sequencing on 131 patients with sporadic
   TOF. By applying frequency-based filtering using both our internal
   control cohort of 6004 individuals without CHDs ([79]16) and public
   databases including Exome Aggregation Consortium (ExAC) and Genome
   Aggregation Database (gnomAD), we defined rare loss-of-function (LoF)
   variants as potential pathogenic candidates and found that these
   variants were significantly enriched in the cilia-related pathways in
   patients with TOF [false discovery rate (FDR)–adjusted P <0.05].
   Further gene-level association analyses revealed that ciliary genes
   exhibited a notably higher burden of rare LoF variants in TOF cases
   compared to controls. The odds ratios (ORs) of ciliary gene variants
   for patients with TOF were calculated with different numbers of LoF
   variants and detrimental missense variants, and the results implied the
   potential oligogenic or multigenic inheritance pattern. Moreover, mouse
   models carrying heterozygous or homozygous knockout/mutant alleles of
   these ciliary genes reproduced a spectrum of cardiac phenotypes
   resembling TOF. Double-heterozygous mice generated through pairwise
   crosses of ciliary gene mutant lines exhibited an even higher incidence
   of TOF-like phenotypes, compared to single-heterozygous or homozygous
   mutants, suggesting additive or synergistic effects between distinct
   ciliary gene variants. In conclusion, this study demonstrates the
   pathogenic role of ciliary gene variants in OFT development and
   provides compelling genetic evidence for the etiology of nonsyndromic
   TOF, with both human and mouse model data supporting a multigenic
   inheritance mechanism underlying its pathogenesis.

RESULTS

Identification of pathogenic variants in known TOF-causative genes

   We first examined rare and deleterious variants in known causative TOF
   genes. To compile a list of candidate genes, we searched the OMIM and
   ClinVar databases and integrated pathogenic genes identified in several
   large-scale, multicenter studies on TOF using whole-exome sequencing
   (WES) and whole-genome sequencing (WGS) ([80]3, [81]10, [82]11, [83]13,
   [84]17, [85]18). This resulted in the inclusion of 11 genes: TBX1,
   NKX2-5, GATA4, ZFPM2, GATA6, JAG1, NOTCH1, FLT4, FOXC2, TBX5, and KDR
   for pathogenicity evaluation.

   We then assessed rare and deleterious variants in these known
   TOF-causative genes (fig. S1A). Rare variants were defined as those
   with a minor allele frequency (MAF) < 0.5% in our study cohorts
   and < 0.5% in both the ALL and East Asian populations of gnomAD and
   ExAC; variant pathogenicity was evaluated through manual review
   following the American College of Medical Genetics and Genomics (ACMG)
   guidelines, or annotation in the ClinVar database (table S1).
   Ultimately, we identified 15 pathogenic and likely pathogenic (P/LP)
   variants across nine known TOF genes, including 2 (13.3%) LoF variants,
   12 (80.0%) missense variants, and 1 (6.7%) start lost variant (fig. S1,
   B to D). These P/LP variants were detected in 16 patients, yielding a
   12.21% contribution to the overall TOF incidence.

Identification of ciliary gene variants in nonsyndromic TOF patients

   Because of the relatively low explanatory power of known TOF-causative
   genes, we then conducted additional case-control analyses at both the
   gene set (pathway) and gene levels to identify enriched genes and
   potential genetic defects associated with nonsyndromic TOF. Since the
   accumulation of mild-effect variants may collectively increase disease
   susceptibility, gene set level analysis offers valuable insight into
   disease mechanisms and may aid in potential gene discovery. Following
   the screening procedure shown in [86]Fig. 1A, we analyzed rare and
   deleterious variants in 131 patients with nonsyndromic TOF. Variants
   were filtered on the basis of MAF_all <0.5% in our TOF and control
   individuals, as well as <0.5% in both of the populations (ALL and East
   Asians population) from gnomAD and ExAC. Deleterious variants were then
   defined as LoF variants, namely frameshift, stop gain, and splice site
   variants (splice acceptor/donor). Last, a total of 2018 rare LoF
   variants on 1777 genes were identified in 131 patients with TOF.
   Subsequently, we performed pathway enrichment analysis and found
   significant enrichment in the cilia-related pathways (P = 3.63 × 10^−8,
   FDR<0.05). A similar result was observed for high-confidence LoF
   variants predicted by the LoF transcript effect estimator ([87]19)
   package ([88]Fig. 1B). There were 108 candidate genes involved in the
   cilium pathway (table S2) with different localizations in cilia (fig.
   S2). In total, 87 of 131 cases harbored ciliary gene LoF variants,
   accounting for 66.4% of TOF cases.

Fig. 1. Identification of ciliary gene variants in patients with nonsyndromic
TOF.

   [89]Fig. 1.
   [90]Open in a new tab

   (A) Flowchart depicting the criteria and results for variant selection
   after next-generation sequencing. Rare and putatively damaging variants
   were filtered according to the MAFs in our reference cohort and public
   databases (ExAC and gnomAD), along with the predicted deleterious LoF
   variants. (B) Genes with rare LoF variants or variants evaluated as
   high-confidence (HC) using the LoF transcript effect estimator (LOFTEE)
   package were both significantly enriched in cilia-related pathways
   (FDR < 0.05). (C) Q-Q plot comparing observed versus expected P values
   for LoF variants in each gene between cases (n = 131) and controls
   (n = 6004), using two-sided Fisher’s exact test. Thirteen genes with
   FDR < 0.2 are highlighted, with ciliary genes marked in red. (D)
   Pathway enrichment analysis of significantly mutated genes that met the
   threshold criteria. (E to I) Immunofluorescence staining of ciliary
   axonemes (green, acetylated α-tubulin), basal bodies (red, γ-tubulin),
   and the nucleus (blue) showing the expression of cilia throughout OFT
   development from E10.5 to E18.5. The white arrows indicate cilia. Three
   levels of magnification were shown from left to right. Scale bars, 200,
   20, and 10 μm. (J) Tomographic electron microscopy of the OFT at E13.5
   showing a straight-cut axoneme, a primary cilium (typical 9+0
   microtubule configuration), and a motile cilium (typical 9+2
   microtubule configuration). Scale bars, 200, 100, 100 nm. (K) Spatial
   expression of all ciliary genes in cilium pathway across different
   cardiac structures. The centerline on the boxplot represents the
   median, the box limits indicate the 25th and the 75th percentiles, and
   the whiskers extend to the minimum and maximum values. The statistical
   analysis was carried out using Wilcoxon–Mann-Whitney test. *P < 0.05,
   **P < 0.01.

   To further explore the candidate pathogenic genes in TOF, we also
   conducted gene-based association analyses in rare deleterious variants.
   We aggregated the qualifying alleles into genes and tested for
   differences between TOF cases and controls using Fisher’s exact test.
   For genes with LoF variants, 13 genes showed a significantly higher
   burden of LoF alleles in TOF cases after multiple testing (FDR < 0.2
   threshold; [91]Fig. 1C), with pathways related to ciliary membrane and
   cilium enriched in these genes (FDR < 0.05; [92]Fig. 1D). Notably,
   three ciliary genes, MOSMO, TAS2R46, and TTLL3, were among the
   significantly mutated genes. In particular, MOSMO is a well-established
   CHD risk gene which is active in ciliary membrane ([93]20). Loss of
   MOSMO results in severe developmental defects in mouse embryos,
   including heterotaxy, skeletal abnormalities, and complex CHDs such as
   double outlet right ventricle (DORV), transposition of the great
   arteries (TGA), and persistent truncus arteriosus, which arise due to
   defects in OFT development ([94]21). In parallel, we performed
   gene-based association tests for predicted detrimental missense
   variants, using multiple in silico tools (AlphaMissense, CADD, SIFT,
   and PolyPhen-2) to evaluate variant pathogenicity. Ciliary genes
   consistently ranked among the top genes with higher burdens of
   detrimental missense variants in TOF cases compared to in-house
   controls across all algorithms (fig. S3). Together, we imply the
   potential role of ciliary gene variants in the pathogenesis of TOF.

Cilia are spatially and temporally regulated during cardiac OFT development

   Since TOF is characterized by abnormal development of the cardiac OFT,
   we initiated experiments to determine whether cilia contribute to OFT
   development by first performing immunofluorescence assays. As shown in
   [95]Fig. 1, E to I, cilia were distributed throughout all stages of
   cardiac OFT development from embryonic day 10.5 (E10.5) to E18.5. These
   findings were confirmed by transmission electron microscopy (TEM),
   which revealed both typical “9+0” (primary cilia) and “9+2” (motile
   cilia) microtubule configurations in the cardiac OFT ([96]Fig. 1J).

   In addition to the OFT, primary cilia were also observed in other
   cardiac regions, such as the atrioventricular canal (AVC), atria,
   ventricles, and valves, at various stages of development (fig. S4, A to
   M). Notably, during early heart development, primary cilia were most
   abundant in the OFT, where they were likely involved in the regulation
   of key signaling pathways, blood flow dynamics, and guiding neural
   crest cell migration, contributing to OFT development and the proper
   separation of the aorta and PA. In contrast, motile cilia exhibited a
   distinct spatiotemporal expression pattern and were less abundant
   compared to primary cilia, with a role primarily in biomechanical
   regulation, contributing to fluid dynamics, extracellular matrix (ECM)
   remodeling, and cell migration, particularly in the OFT and valve
   regions (fig. S4, N to Z).

   To further investigate the heterogeneity in the distribution of cilia
   within distinct anatomical regions of the heart, by intersection with
   spatial transcriptomics data of E10.5 hearts ([97]22), we revealed that
   ciliary genes showed notably higher spatial expression in the OFT
   [including in OFT cardiac neural crest cells (cNCCs), OFT
   cardiomyocytes, and mesenchymal cells] (P < 0.05), supporting the
   embryological hypotheses of aberrant positioning of the OFT in TOF
   ([98]Fig. 1K). In parallel, we performed single-cell RNA sequencing
   (scRNA) on anatomically dissected E10.5 hearts, dividing them into six
   regions: OFT, AVC, left atrium (LA), right atrium (RA), left ventricle
   (LV), and right ventricle (RV). Consistent with the spatial
   transcriptomics results, ciliary gene expression was highest in the AVC
   and OFT regions (fig. S5, A to B). Given the role of the AVC in
   endothelial-to-mesenchymal transition and cushion formation, which
   gives rise to both AV and OFT cushions, this expression pattern further
   highlights the regional specificity of ciliary gene activity during OFT
   development.

   To provide broader context for ciliary gene expression, we also
   analyzed transcriptomic datasets across multiple tissues ([99]23) and
   cardiac cell types. Bulk RNA-seq data showed that ciliary genes are
   broadly expressed during development, with notably high expression in
   the heart, brain, and testes. Genes highly expressed in the heart
   (clusters 5, 6, and 7) exhibited a greater LoF variant burden in TOF
   cases compared to controls (OR, 1.96, P = 0.002, 95% CI, 1.255 to
   2.976), supporting their potential contribution to cardiac-specific
   phenotypes (fig. S5C). ScRNA-seq analysis further revealed that ciliary
   genes are expressed across diverse cardiac cell types, including
   cardiomyocytes, mesenchymal cells, and endocardial cells, where they
   participate in distinct signaling pathways critical for OFT development
   (fig. S5D). Together, these findings highlight that although ciliary
   genes are broadly expressed across multiple tissues and cardiac
   regions, they show particularly elevated expression and functional
   relevance in the OFT during early heart development, supporting their
   critical role in the pathogenesis of TOF.

Evaluation of the impacts of ciliary gene variants in the pathogenesis of
nonsyndromic TOF

   Increasing evidence indicates that in addition to monogenic
   inheritance, CHD may involve oligogenic or multigenic inheritance,
   which is characterized by the presence of damaging heterozygous
   mutations across more than one gene ([100]24, [101]25). To
   systematically evaluate the overall contribution of ciliary gene
   variants to TOF pathogenesis and assess their predominant inheritance
   pattern, we performed gene set–based association analyses focusing on
   deleterious variants in all genes involved in the cilium pathway
   ([102]Fig. 2A and table S3). ORs were calculated by comparing the
   number of patients with TOF and controls carrying different burdens of
   LoF variants within this predefined gene set. We observed an increased
   burden of ciliary gene LoF variants among patients with TOF compared to
   controls (OR, 1.540; 95% CI, 1.056 to 2.775; P = 0.0205). This
   association was more pronounced in individuals carrying two or more LoF
   variants (OR, 1.888; 95% CI: 1.265 to 2.778; P = 0.0017; [103]Fig. 2B).
   A similar enrichment pattern was observed for putatively detrimental
   missense variants in patients with TOF, with an OR of 1.590 (95% CI:
   1.106 to 2.284; P = 0.0105) in individuals carrying more than two such
   variants ([104]Fig. 2B). When combining the impacts of LoF variants and
   detrimental missense variants, we observed that as the number of
   variants increased, the risk for TOF increased accordingly.
   Specifically, when the number of variants reached six, our TOF cohort
   exhibited greater than threefold enrichment (OR, 3.158; 95% CI, 1.381
   to 6.371; P = 0.0038), thus supporting a multigenic etiology of ciliary
   genes in the pathogenesis of TOF ([105]Fig. 2B).

Fig. 2. Rare LoF variants and detrimental missense variants burden analyses
revealed the multigenic inheritance effect of ciliary genes in patients with
nonsyndromic TOF.

   [106]Fig. 2.
   [107]Open in a new tab

   (A) Flowchart illustrating the selection of ciliary gene variants for
   gene set–based association analyses between 131 patients with
   nonsyndromic TOF and 6004 in-house control individuals. (B) Forest plot
   showing the burden test results for patients with TOF carrying
   different thresholds of LoF variants, detrimental missense variants,
   and combined deleterious variants (LoF plus detrimental missense
   variants) in ciliary genes, compared to control individuals. The number
   of participants in each category is displayed, where N represents the
   number of distinct heterozygous variants detected across different
   ciliary genes. The analysis includes biallelic heterozygous variants.

   To control for potential confounding effects from biallelic
   inheritance, we reevaluated the model after excluding individuals
   carrying biallelic variants. Among the 131 patients with TOF, such
   variants were identified in 6.87% (9 of 131), with LoF variants
   accounting for 2.29% (3 of 131) and detrimental missense variants for
   3.05% (4 of 131) (fig. S6A). Notably, a substantial proportion (up to
   66.41%) of TOF cases carried more than two pathogenic variants. After
   removing biallelic individuals, the multiheterozygous burden analysis
   in the remaining 122 TOF cases and 6004 controls yielded consistent
   trends, with increasing TOF risk observed as the number of pathogenic
   variants increased. Specifically, ORs were 2.136, 2.601, and 3.149 for
   individuals carrying more than four, five, and six variants,
   respectively (fig. S6B). These results reinforce a multigenic
   inheritance model in TOF and suggest that cumulative variant burden in
   ciliary genes may contribute to a substantial proportion of unsolved
   cases.

Loss of ciliary genes impairs normal cardiac OFT development

   To explore the causative effects of ciliary genes on cardiac OFT
   development in vivo, we prioritized genes harboring rare LoF variants
   enriched in the cilium pathway ([108]Fig. 1A and fig. S7A). Among 90
   previously unreported candidate genes (i.e., genes not confirmed to be
   associated with CHD in OMIM database or genes without CHD phenotypes in
   mice from MGI database), genes with recurrent variants were selected
   (n = 17), and after filtering out variants that were poorly conserved
   (<60%) in humans and mice (n = 6), failed Sanger sequencing validation
   (n = 1), and had technical difficulties in mouse model construction
   (n = 4), we ultimately opted to construct mouse models for the retained
   six ciliary genes, including AGBL2, MLF1, PPEF2, PKHD1L1, SPTBN5, and
   TEKT3. Together, we identified a heterozygous frameshift variant
   (c.2683_2684delTC, p.Ser895fs) and a heterozygous stop-gain variant
   (c.742C>T, p.Arg248*) in AGBL2; a heterozygous stop-gain variant
   (c.778C>T and p.Arg260*) in two individuals in MLF1; three heterozygous
   frameshift variants (c.1820_1821delTT, p.Phe607fs; c.2714_2717dupCAGA,
   p.Glu906fs; and c.4727delG, p.Gly1576fs) in PKHD1L1; a heterozygous
   stop-gain variant (c.1038G>A, p.Trp346*) and a heterozygous frameshift
   variant (c.1632_1633dupAC, p.Leu545fs) in PPEF2; three heterozygous
   frameshift variants (c.10970_10973delATGA, p.Asn3657fs; c.9417dupT,
   p.Val3140fs; c.1318_1319delAG, p.Ser440fs) in SPTBN5; a heterozygous
   stop-gain variant (c.544G>T, p.Glu182*) and a heterozygous frameshift
   variant (c.546_547delGCinsT, p.Glu182fs) in three individuals in TEKT3
   (fig. S7B and table S4).

   Using CRISPR-Cas9 technology, we then generated six ciliary gene
   knockout or mutant mouse models (fig. S8). Sanger sequencing of
   homozygous mutant mice confirmed the presence of deletions in Agbl2,
   Pkhd1l1, Sptbn5, Tekt3, and anthropomorphic variants in Mlf1 and Ppef2
   (table S5), which were predicted to cause premature translational
   termination or gene inactivation. Meanwhile, no predicted off-target
   sites (OTSs) were found (fig. S9 and table S6). Compared with wild-type
   mouse hearts, protein immunoblotting and immunofluorescence showed
   substantial decrease in the expression of the corresponding proteins in
   the heart tissues of knockout or mutant mice ([109]Fig. 3A and figs.
   S10 to S14A). As homozygous inactivation of mouse ciliary genes has
   been reported to cause embryonic lethality ([110]26, [111]27),
   postnatal day one (P1) mice from heterozygous parents were examined.
   The genotype ratio from P1 mice was in accordance with the expected
   Mendelian ratio among the offspring ([112]Fig. 3B and figs. S10 to
   S14B), whereas homozygous P1 mice suffered from higher postnatal
   lethality ([113]Fig. 3, C to D, and figs. S10 to S14, C and D). Then,
   we examined the cardiac phenotypes of each strain. Dissected and fixed
   P1 hearts showed abnormal and malrotation of the aorta and PA in
   Pkhd1l1^−/− hearts ([114]Fig. 3E). Further experiment of hematoxylin
   and eosin (H&E) staining of P1 Pkhd1l1^−/− hearts revealed cardiac
   malformations in a subset of homozygous mice. Among them, 9 of 73
   (12.33%) exhibited an OA in combination with a VSD (OA and VSD), while
   7 of 73 (9.59%) showed muscular VSDs (MU-VSD), and 11 of 73 (15.07%)
   presented with secundum atrial septal defects (ASD), together
   accounting for 39.7% of the homozygotes. Notably, cardiac defects were
   also observed in heterozygous mice, including OA and VSD in 6 of 94
   (6.38%), MU-VSD in 5 of 94 (5.32%), ASD in 5 of 94 (5.32%), and
   combined MU-VSD and ASD in 1 of 94 (1.06%) heterozygotes. However, no
   similar structural defects were observed in P1 hearts of wild-type mice
   ([115]Fig. 3, C and F). Mice that had died on P1 day mostly exhibited
   structural cardiac defects ([116]Fig. 3D). Moreover, the PA widths of
   the hearts in P1 Pkhd1l1^−/− mice appeared to be thinner than those of
   the normal controls ([117]Fig. 3G). Together, knockout of the ciliary
   gene Pkhd1l1 resulted in a spectrum of cardiac malformations, including
   OA with VSD and other defects affecting septal development and OFT
   formation, particularly VSD.

Fig. 3. Pkhd1l1-knockout mice showed typical abnormal cardiac OFT phenotypes.

   [118]Fig. 3.
   [119]Open in a new tab

   (A) Immunofluorescence staining of PKHD1L1 in E10.5 OFT sections from
   wild-type (Pkhd1l1^+/+) and homozygous (Pkhd1l1^−/−) mice. Scale bars,
   200 and 10 μm. (B) Genotype frequency analysis of offspring from
   Pkhd1l1^+/− heterozygote crosses, showing the expected Mendelian ratio
   for P1 newborns (χ2 = 0.2818, P>0.05). (C) Postnatal lethality and
   cardiac defects (OA and VSD, MU-VSD, ASD, and MU-VSD with ASD) in
   wild-type, heterozygous, and homozygous P1 offspring. (D) Proportions
   of normal and defective phenotypes in Pkhd1l1 knockout mice that died
   on P1. (E) Malrotation of the aorta and PA in Pkhd1l1^−/− hearts at P1.
   Scale bar, 1000 μm. (F) Representative H&E staining of abnormal hearts
   from Pkhd1l1^−/− homozygous offspring at P1, showing OA and VSD,
   MU-VSD, and ASD (asterisk). Scale bar, 500 μm. (G) Measurements of P1
   hearts indicated thinner pulmonary arteries due to Pkhd1l1 deletion.
   The outline and diameter of the PA were indicated by black lines.
   Statistical analysis was performed using an unpaired two-tailed
   Student’s t test (**P < 0.01). n>30 mice per group. Scale bar, 1000 μm.
   (H) E10.5 wild-type and Pkhd1l1 deletion embryos were shown, with
   enlarged view of the OFT below. Homozygous mutants exhibited reduced
   OFT length compared to littermate controls, confirmed by quantitative
   measurements normalized by dorsal-ventral axis lengths. The OFT length
   and dorsal-ventral axes length are represented by white lines and red
   lines, respectively. Statistical analysis was performed using an
   unpaired two-tailed Student’s t test (*P < 0.05). n >15 mice per group.
   Scale bar, 1000 μm. (I and J) Tomographic electron microscopy of E13.5
   cardiac OFT cross sections revealed typical 9+0 microtubule structures
   in Pkhd1l1^+/+ mice, while Pkhd1l1^−/− mice exhibited severely
   disorganized 9+0 axonemes. Scale bars, 200 nm. OA, overriding aorta;
   MU-VSD, muscular ventricular septal defect; ASD, secundum atrial septal
   defect.

   Proper folding and development of the heart tube requires recruitment
   of second heart field (SHF) cardiac precursors that contribute to the
   lengthening of the heart tube, eventually giving rise to much of the
   right ventricle, atria, and OFT ([120]28). Therefore, to further
   investigate the impact of ciliary gene deletion on early OFT
   development and morphogenesis, we examined E10.5 embryos and observed
   decreased size and inadequate looping of the OFTs in several homozygous
   Pkhd1l1^−/− mice ([121]Fig. 3H). To account for potential variations in
   overall embryonic development, we normalized OFT lengths to
   dorsal-ventral axis lengths. Blinded quantitative analysis revealed a
   significant reduction in normalized OFT length (P < 0.05) in
   Pkhd1l1^−/− embryos compared to littermate controls, indicating that
   ciliary gene loss leads to intrinsic defects in OFT elongation during
   early cardiac development ([122]Fig. 3H).

   Consistent with the phenotypes observed in Pkhd1l1^−/− mouse hearts,
   knockout or mutation of other ciliary genes also led to similar cardiac
   phenotypes, with a phenotypic incidence of approximately 30% across
   strains. In Agbl2^−/− mice, 31.9% (15 of 47) of homozygotes exhibited
   cardiac defects, including OA and VSD in one mouse (2.13%), MU-VSD in
   eight (17.02%), ASD in five (10.64%), and combined MU-VSD and ASD in
   one (2.13%). In Mlf1^KI/KI mice, the defect rate was 23.4% (15 of 64),
   comprising OA and VSD in one (1.56%), MU-VSD in four (6.25%), ASD in
   nine (14.06%), and MU-VSD with ASD in one (1.56%). In Ppef2^KI/KI mice,
   28.1% (25 of 89) displayed cardiac defects, including OA and VSD in 5
   (5.62%), MU-VSD in 5 (5.62%), and ASD in 15 (16.85%). In Sptbn5^−/−
   mice, 23.8% (15 of 63) had malformations, with OA and VSD in five
   (7.94%), MU-VSD in three (4.76%), ASD in six (9.52%), and combined
   MU-VSD and ASD in one (1.59%). In Tekt3^−/− mice, 31.3% (15 of 48) were
   affected, including OA and VSD in one (2.08%), MU-VSD in five (10.42%),
   ASD in eight (16.67%), and MU-VSD with ASD in one (2.08%) (figs. S10 to
   S14, C and D). In addition, we identified OFT and septal malformations
   in each heterozygous strain, whereas no cardiac developmental defects
   were observed at P1 in wild-type mouse hearts in ciliary genes above,
   which illuminated the supposed burden and pathogenicity of patients
   with heterozygous variants (figs. S10 to S14E). Similar results were
   observed in comparisons of PA widths and E10.5 normalized OFT lengths,
   with both parameters significantly reduced in all knockout and mutant
   mouse hearts (P < 0.05) (figs. S10 to S14, F and G). Overall, these
   findings implicate an important role of ciliary gene variants in OFT
   development.

   To assess whether mice with combined OA and VSD exhibit a TOF-like
   phenotype, we further measured PA widths in these P1 mice. Given the
   lack of an established reference range for PA widths in P1 mice, we
   compiled measurements from wild-type P1 littermates derived from
   heterozygous crosses of each ciliary gene knockout/mutant strain,
   resulting in an average PA width of 415.23 μm, which was used as the
   normal reference for P1 mice. PA widths in both heterozygous and
   homozygous ciliary gene knockout/mutant mice with combined OA and VSD
   were then compared against this reference. All affected mice exhibited
   narrower pulmonary arteries than the wild-type average, with a mean PA
   width of 371.52 μm (fig. S15). The simultaneous presence of OA, VSD,
   and PA narrowing in these P1 mice recapitulates the key features of the
   human TOF phenotype. Therefore, the incidence of TOF-like defects
   across different strains is as follows: Agbl2^−/−: 2.13%, Mlf1^KI/KI:
   1.56%, Pkhd1l1^−/−: 12.33%, Ppef2^KI/KI: 5.62%, Sptbn5^−/−: 7.94%, and
   Tekt3^−/−: 2.08%.

   Building on the observed phenotypes in single ciliary gene
   knockout/mutant mice, we next sought to investigate the potential
   cumulative effects of compound heterozygosity. We performed pairwise
   crosses of our ciliary gene knockout/mutant mice, generating five
   distinct double-heterozygous lines ([123]Fig. 4A). A total of 122 P1
   hearts were collected for cardiac defect phenotype analysis. The
   overall incidence of cardiac defects in the double-heterozygous mice
   was 33.6%, with the incidences for each strain as follows:
   Pkhd1l1+Tekt3 at 17.4%, Tekt3+Ppef2 at 50.0%, Ppef2+Mlf1 at 42.1%,
   Pkhd1l1+Mlf1 at 47.6%, and Tekt3+Sptbn5 at 28.6% ([124]Fig. 4A). Among
   these, OA was observed in 11.48% (14 of 122) of cases, while the
   incidences of perimembranous VSD (PM-VSD), MU-VSD, and ASD were 4.10,
   11.48, and 8.20%, respectively ([125]Fig. 4B). The increased incidence
   of OA in double-heterozygous mice further suggests a potential
   cumulative or interactive genetic effect. Moreover, 21 hearts (17.21%)
   were collected from mice that had died on P1 day, with a markedly
   higher defect rate of 61.90%. In this group, 19.05% (4 of 21) exhibited
   OA, 9.52% had PM-VSD, 28.57% had MU-VSD, and 14.29% had ASD ([126]Fig.
   4C), indicating that early postnatal mortality was closely associated
   with cardiac abnormalities. A comparative analysis of overall cardiac
   defect incidence revealed consistently higher frequencies in
   double-heterozygous mice compared to single-heterozygous strains
   ([127]Fig. 4D). In addition, the incidence of OA was higher in
   double-heterozygous mice compared to single-heterozygous strains and,
   in most cases, exceeded that observed in homozygous mice ([128]Fig. 4,
   E to F). In conclusion, the findings from both our human cohort and
   mouse models support a multigenic inheritance model in TOF. The
   incomplete penetrance of cardiac phenotypes in heterozygous mice
   further suggests that the pathogenesis of TOF may arise from the
   combined effects of deleterious heterozygous variants on more than one
   ciliary gene.

Fig. 4. Cardiac phenotypic analysis of double-heterozygous ciliary gene
knockout/mutant mice.

   [129]Fig. 4.
   [130]Open in a new tab

   (A) Breeding scheme for generating double-heterozygous mice. The
   numbers in colored boxes below indicate the number of cardiac defects
   observed in each strain, with each color representing a specific
   phenotype: red for OA + VSD, purple for OA + VSD + ASD, green for
   PM-VSD, blue for MU-VSD, and orange for ASD. (B) The number and
   proportion of each cardiac defect phenotype in the double-heterozygous
   mice. (C) The number and proportion of each cardiac defect phenotype in
   the double-heterozygous mice that died before heart collection. (D)
   Comparison of overall cardiac defect probability between
   single-heterozygous and double-heterozygous ciliary gene
   knockout/mutant mice. (E) Comparison of OA probability between
   single-heterozygous and double-heterozygous ciliary gene
   knockout/mutant mice. (F) Representative H&E staining of abnormal
   hearts from double-heterozygous mice at P1, showing phenotypes of OA,
   PM-VSD, MU-VSD, and ASD (asterisk). Scale bars, 500, 500, and 1000 μm.
   PM-VSD, perimembranous VSD.

Loss of ciliary genes results in aberrant ciliogenesis during cardiac OFT
development

   Given the well-established role of cilia structure in maintaining
   normal ciliary function ([131]29), we further investigated the effects
   of variants in ciliary genes on cilia phenotype and function in OFTs.
   We first examined the subcellular localization of these ciliary genes
   in both mouse embryonic fibroblasts (MEFs) and E10.5 mouse OFTs.
   Immunofluorescence staining analysis revealed that Mlf1 and Pkhd1l1
   colocalized with acetylated α-tubulin and specifically accumulated in
   cilia axonemes; Sptbn5 colocalized with both acetylated α-tubulin and
   γ-tubulin along the entire cilia; Ppef2 and Tekt3 were localized to
   basal bodies of cilia; and Agbl2 was located at the base of cilia and
   overlapped with Nek2, a marker of centriole ([132]Fig. 5A and fig.
   S16).

Fig. 5. Loss of ciliary genes disrupted ciliogenesis during cardiac OFT
development.

   [133]Fig. 5.
   [134]Open in a new tab

   (A) Schematic of a cilium showing structural features and selected
   ciliary genes. Ciliary genes (red) were immunostained with axonemes
   (acetylated α-tubulin), basal bodies (γ-tubulin), and centrioles
   (NEK2). Nuclei were stained with DAPI (blue). Mlf1 and Pkhd1l1
   localized to cilia axonemes; Sptbn5 localized along the entire cilia;
   Ppef2 and Tekt3 localized to basal bodies of cilia; and Agbl2 was
   observed both at the basal body and centriole of cilia. Boxed and
   color-coded numbers following each gene indicate the number of cardiac
   defects observed in knockout/mutant mouse hearts: red for OA + VSD,
   purple for OA + VSD + ASD, blue for MU-VSD, orange for ASD, and green
   for MU-VSD + ASD. Scale bars, 20 μm. (B) Immunofluorescence staining of
   cilia in the cardiac OFTs of ciliary gene knockout/mutant mice and
   control mice at E13.5. Axonemes (acetylated α-tubulin, green) and basal
   bodies (γ-tubulin, red) were shown. The right column showed higher
   magnification of the boxed area from the left. Scale bars, 20 and 10
   μm. (C) Box plots of cilia prevalence (left) and cilia length (right)
   in the cardiac OFTs of ciliary gene knockout/mutant mice and control
   mice at E13.5. (Left): Error bars represent means ± SEM. (Right): The
   blue boxes showed the lower, median, and upper quartiles of
   measurements. Each black dot represented the cilia length from a single
   cell. Homozygous knockout/mutant mice exhibited shorter cilia and a
   lower cilia percentage compared to wild-type control mice, indicating
   disrupted ciliogenesis. Statistical significance was assessed by an
   unpaired two-tailed Student’s t test or Wilcoxon–Mann-Whitney test when
   appropriate. *P < 0.05, **P < 0.01, ***P < 0.001. n = 10 to 20 OFT
   sections from five hearts for each group were measured percentage of
   cilia and cilia length across all OFT sections (n > 400 cilia per
   group).

   Correct regulation of cilia length is essential for proper cilium
   biogenesis and function ([135]30), whereas defects in ciliary axonemes
   or basal proteins may lead to abnormal cilia length ([136]31). To
   assess the impact of ciliary gene variants on ciliogenesis during OFT
   development, we collected E13.5 hearts from wild-type and
   knockout/mutant mice. Immunofluorescence examination of E13.5 OFTs
   revealed significantly reduced cilia number and shorter lengths in
   homozygous mice compared with wild-type littermates (P < 0.05)
   ([137]Fig. 5, B and C). Ultrastructural analysis using TEM further
   revealed that Pkhd1l1^−/− mouse OFT tissues displayed severely
   disorganized ciliary axonemes, with microtubule loss and irregular
   ciliary distribution and orientation, compared to the typical 9+0
   microtubule arrangement in wild-type Pkhd1l1^+/+ mice ([138]Fig. 3, I
   and J). Thus, these findings highlight the critical role of ciliary
   genes in ciliogenesis and their essential function located in cardiac
   OFT development.

Cilia orchestrate OFT development through dysregulation of sonic Hh signaling

   Cilia serve as a transduction center for a series of
   development-related signaling pathways during cardiovascular
   development, including Hedgehog (Hh) signaling pathway, WNT signaling
   pathway, and transforming growth factor–β (TGF-β)/bone morphogenetic
   protein (BMP) signaling pathway ([139]32). Both abnormalities in cilia
   structure and changes in the number of cilia may cause alterations in
   signaling pathways, which in turn lead to a range of developmental
   disorders. To further explore the shared mechanisms underlying OFT
   abnormalities caused by ciliary gene knockout or mutation, we carried
   out whole-transcriptome sequencing (RNA-seq) of E10.5 hearts from
   wild-type and knockout/mutant littermates of each strain. Pathway
   enrichment analysis from Kyoto Encyclopedia of Genes and Genomes (KEGG)
   database indicated differentially expressed genes (DEGs) were
   associated with cellular communication, including cell adhesion
   molecules, tight junction, and ECM-receptor interaction, as well as
   pathways related to cardiac muscle contraction and Hh signaling.
   Intriguingly, the Hh signaling pathway was the only pathway
   consistently and significantly enriched across all ciliary gene
   knockout and mutant strains (FDR < 0.05) ([140]Fig. 6A and fig. S17),
   providing strong support for its critical role in mediating early
   cilia-dependent cardiac development. To further investigate the
   specific cell type responsible for the observed phenotype resulting
   from cilia abnormalities, we then evaluated the expression of cardiac
   marker genes in each group. The results showed reduced expression of
   several myocardial marker genes Tnnt2 and Myl2 in ciliary gene
   knockout/mutant groups ([141]Fig. 6B), whereas markers of other cardiac
   cell remained unchanged (fig. S18A). Meanwhile, as shown in [142]Fig. 6
   (C and D), a reduced number of cilia was also observed in
   cardiomyocytes of knockout/mutant mouse OFTs from same-litter crosses
   of each strain. Then, by integrating scRNA-seq data ([143]Fig. 6E and
   fig. S18, B and C) and applying gene set variation analysis (GSVA), we
   further evaluated Hh pathway activity and found lower activities of Hh
   pathway in cardiomyocytes of ciliary gene knockout/mutant strains
   ([144]Fig. 6F). In addition, DEGs in these cardiomyocytes were enriched
   in pathways in heart development, circulatory system development, and
   cardiac muscle tissue development ([145]Fig. 6G), suggesting the
   abnormal cardiomyocyte function and cardiac development affected by
   cilia abnormalities. The activity of Hh signaling pathway is mainly
   transported downward through the activated membrane protein receptor
   SMO (Smoothened) and then activates the nuclear transcription factor
   GLI family (including GLI1, GLI2, and GLI3), which subsequently enters
   the nucleus and initiates the expression of target genes downstream
   ([146]33). Western blot assays showed that expression of GLI1 and GLI2
   was notably decreased in the hearts of knockout and mutant mice (fig.
   S18D), confirming the disruption of Hh pathway activity.

Fig. 6. Integrated analysis of bulk RNA-seq and scRNA-seq correlated cilia
loss with dysregulation of Hh signaling pathway and decreased Gata4 and
Nkx2-5 expression.

   [147]Fig. 6.
   [148]Open in a new tab

   (A) The Sankey diagram showing the top 10 enriched pathways in the
   hearts of homozygous ciliary gene knockout/mutant mice versus
   littermate controls at E10.5. (B) The expression of myocardial marker
   genes in E10.5 hearts of wild-type and ciliary gene knockout/mutant
   groups. Bars represent means ± SEM, with statistical analysis by
   unpaired two-tailed Student’s t test or Wilcoxon–Mann-Whitney test.
   *P < 0.05, **P < 0.01, ***P < 0.001. n = 3 to 4 biological replicates
   per group, with 10 hearts per replicate. (C) Immunofluorescence
   staining for cilia in cardiomyocytes of wild-type and knockout/mutant
   mouse OFTs from the same-litter mice, with cardiac troponin t (cTnT) in
   green and axonemes (acetylated α-tubulin) in red. Scale bars, 200, 20,
   and 10 μm. n = 3 to 4 replicates per group. (D) Quantification of cilia
   percentage in cardiomyocytes in (C), presented as means ± SEM.
   Statistical significance was assessed using unpaired two-tailed
   Student’s t test or Wilcoxon–Mann-Whitney test. ***P < 0.001. (E)
   Uniform manifold approximation and projection (UMAP) clustering of
   cells, showing distinct clusters. (F) The Hh pathway activity in the
   cardiomyocyte cluster of wild-type and each ciliary gene
   knockout/mutant group. Statistical analysis was performed using
   unpaired two-tailed Student’s t test or Wilcoxon–Mann-Whitney test.
   *P < 0.05, **P < 0.01, ***P < 0.001. (G) GO analyses of the scRNA-seq
   datasets with DEGs (FDR < 0.05) in cardiomyocytes of knockout/mutant
   strains versus wild-type OFTs. The top 12 enrichment pathways shown in
   bubble charts, with cardiac development pathways marked in red. (H)
   Flow chart for selecting cardiac target genes regulated by
   cilia-mediated Hh signaling. (I) The down-regulated genes in all cilia
   knockout/mutant strains expressed by SHF and/or cNCCs. (J) IGV tracks
   showing Gli1 and Gli2 locations at Gata4 and Nkx2-5 locus from CUT&Tag
   analysis in E10.5 mouse hearts.

   The Hh signaling pathway is critical for cardiac OFT development and
   septation. Several studies have reported that Hh signaling pathway is
   required for the survival and migration of cNCCs to the OFT cushions
   and plays an essential role in SHF-mediated OFT elongation and
   septation after cushion formation ([149]34, [150]35). On the basis of
   these findings, we hypothesize that ciliary gene mutations impair Hh
   signaling, affecting gene expression of cardiac progenitor cells in SHF
   and/or cNCCs, thus contributing to abnormal OFT development. Through
   literature review, we identified several genes co–down-regulated in all
   cilia knockout and mutant strains, including SHF-derived genes, namely
   Foxc2, Gata4, Hopx, Nkx2-5, and Srf ([151]36–[152]38), as well as
   cNCC-expressed genes like Cdh2, Foxc2, and Vangl1 ([153]39). Besides,
   we examined the Gli-target genes in heart through performing cleavage
   under targets and tagmentation (CUT&Tag) experiments on E10.5 mouse
   hearts using the GLI1 and GLI2 antibodies. Between them, we found two
   overlapping genes, Gata4 and Nkx2-5 ([154]Fig. 6, H and I), which have
   also been demonstrated to play important roles in cardiac OFT
   morphogenesis and development ([155]40, [156]41). Both Gli1 and Gli2
   showed notable enriched at Gata4 promoter region (≤1 kb); Gli1 binding
   was also observed at promoter (≤1 kb) of Nkx2-5; while most Gli2
   binding sites were around the intergenic region of Nkx2-5 ([157]Fig.
   6J). Motif analysis of Gli1 and Gli2 revealed overrepresentation of
   transcription factor binding sites in Gata4 and Nkx2-5 regions
   respectively (P < 0.05) (fig. S18E). Moreover, immunoblotting analysis
   for GATA4 and NKX2-5 revealed notably reduced protein levels in ciliary
   gene knockout and mutant hearts (fig. S18F). Importantly, a deficiency
   in SHF-derived cardiomyocytes addition leads to a shorter OFT, which is
   consistent with the phenotype we observed in E10.5 mice with ciliary
   gene knockout/mutation ([158]Fig. 3H and figs. S10 to S14G). Together,
   our data indicate that ciliary genes knockout/mutation may result in
   cilia abnormalities in cardiomyocytes, down-regulating the activities
   of Hh signaling pathway and the expression of the Gli-target genes
   Gata4 and Nkx2-5, which ultimately contribute to the abnormal OFT
   development and TOF phenotype.

DISCUSSION

   In this study, we provide a detailed characterization of genetic
   landscape for nonsyndromic TOF. For known causative TOF genes, we
   adopted Clinvar and ACMG standards to classify variant pathogenicity
   and identified 15 P/LP variants in these genes, yielding a 12.21%
   contribution to TOF incidence (16 of 131 cases). Subsequent association
   analyses comparing TOF cases with our in-house controls (6004
   individuals without TOF) identified cilia-related pathways and ciliary
   genes enriched with rare LoF variants both at gene set (pathway) level
   and gene level. Multihit ciliary genes were notably enriched among
   patients with TOF with LoF variants and putatively detrimental missense
   variants, with an OR greater than threefold when the number of combined
   deleterious variants reached six. Through construction of
   knockout/mutant mouse models, including double-heterozygous models, we
   further recapitulated the TOF-like phenotypes observed in patients with
   TOF and explored the potential molecular mechanism. The LoF mouse
   models showed TOF-like phenotypes, a severe reduction in the number and
   length of cilia, and disorganized axoneme arrangements in the cardiac
   OFTs. The incidence of combined OA and VSD was substantially higher in
   double-heterozygous mice, suggesting a cumulative genetic effect. By
   integrating RNA-seq and single-cell sequencing data from cardiac OFTs,
   followed by CUT&Tag analysis in hearts, we speculated that ciliary
   genes knockout/mutation might down-regulate Gli expression by
   disrupting Hh signaling pathway mainly in cardiomyocytes. The target
   genes of Gli (Gata4 and Nkx2-5) in the heart were experimentally
   verified, which in turn lead to reduced OFT length and ultimately TOF
   phenotype ([159]Fig. 7).

Fig. 7. Proposed mechanism of ciliary gene knockout/mutations leading to
tetralogy of Fallot.

   [160]Fig. 7.
   [161]Open in a new tab

   Population-based genetic analyses and mouse models collectively
   highlight the pathogenic role of ciliary genes in cardiac OFT
   development. Knockout or mutation of ciliary genes resulted in reduced
   cilia number and length, as well as disorganized axonemal structure in
   the cardiac OFT. Mechanistically, ciliary dysfunction inactivated Hh
   signaling in cardiomyocytes, leading to transcriptional down-regulation
   of its target genes, Gata4 and Nkx2-5, ultimately disrupting OFT
   morphogenesis and contributing to the development of TOF.

   The ethnic differences between European and Asian populations may lead
   to certain variances in TOF-associated genetic susceptibility. To
   elucidate the discrepancies between our findings and those of the
   European TOF cohort, we examined ultrarare variants in NOTCH1 and FLT4,
   two top TOF-associated mutated genes reported with increased mutational
   burden in TOF cases of European ancestry from prior WES studies
   ([162]13). We then examined the overrepresentation of variants in these
   two genes in our TOF cohort. Among our 131 TOF cases, we detected four
   patients with ultrarare damaging variants: three in NOTCH1 (2.29%) and
   one in FLT4 (0.76%), accounting for 3.05% of the TOF cases. In
   contrast, our control group (n = 6004) harbored 68 variants across 80
   individuals, a rate of 1.33%. Thus, our Chinese patients with TOF had a
   twofold increase in mutational burden compared to controls, although no
   statistically significant association was detected (OR, 2.332; 95%
   CI, 0.611 to 6.354; P = 0.1042). Of note, consistent with findings from
   another Chinese TOF cohort by Tang et al. ([163]42), the proportion of
   carriers in the Chinese population was considerably lower than that of
   the European populations (i.e., 4.46% for NOTCH1, 2.53% for FLT4) (fig.
   S19). Hence, ethnic differences of the population may partly explain
   the diversity of pathogenic mechanisms of TOF.

   Previous studies have reported that cilia are closely associated with
   cardiac development ([164]43). Mouse studies revealed a wide
   distribution of primary cilia throughout the embryonic mouse heart,
   among the ventricular surface of trabeculations in both ventricles and
   were more prominent in the atrial endocardial layer, endocardial
   cushion mesenchyme, and epicardium ([165]44). A large-scale forward
   genetic screening study mutagenized with ethylnitrosourea (ENU)
   constructed a number of mutant mouse models and identified
   approximately 100 variants in 61 genes induced cardiovascular
   phenotypes. This included 34 cilia-related genes, 16 genes involved in
   cilia-transduced cell signaling, and 10 genes regulating vesicular
   trafficking, a pathway important for ciliogenesis and cell signaling
   ([166]32). Although this mutagenesis screen provided compelling
   evidence that cilia may play a central role in CHD pathogenesis, the
   relevance of this finding to human CHD has yet to be fully clarified.
   Moreover, the roles of ciliary gene variants contributing to the
   pathogenesis of TOF have not been studied in a systematic manner,
   either through population-based studies or through mouse models. In the
   present study, we systematically screened rare LoF variants as
   potential pathogenic variants in nonsyndromic TOF. Through gene
   set–based pathway analysis and gene-based association analysis, we
   identified notable enrichment of cilia-related pathways and ciliary
   genes and verified the overall impacts of deleterious ciliary gene
   variants on TOF. To our knowledge, this study represents comprehensive
   investigation of ciliary gene involvement in the pathogenesis of
   nonsyndromic TOF. By applying stringent filtering criteria, we further
   identified six ciliary genes (AGBL2, MLF1, PKHD1L1, PPEF2, SPTBN5, and
   TEKT3) with recurrent rare LoF variants enriched in patients with TOF,
   and conducted corresponding mouse models for subsequent functional
   validation.

   Increasing evidence suggests that a large proportion of CHD might be
   oligogenic or polygenic rather than single genetic inheritance or
   caused by a regulatory network and shared signaling pathways
   constructed by interrelated disease-related genes ([167]45), which
   could account for the reduced penetrance and variable expressivity
   commonly observed in CHD ([168]46). In support of this, an exome-based
   case-control analysis of 249 patients supported a complex and
   heterogeneous genetic architecture of TGA, with 33% of the TGA probands
   most likely affected by oligogenic or polygenic inheritance ([169]24).
   Gifford et al. ([170]47) further demonstrated an oligogenic basis for
   CHD, showing that combined heterozygous mutations in three genes (MKL2,
   MYH7, and NKX2-5) resulted in left ventricular noncompaction in both
   humans and mice. In our study, we observed that heterozygous ciliary
   gene knockout/mutant mice developed mild but discernible cardiac
   phenotypes, with some strains presenting OA and associated defects such
   as ASD and MU-VSD, supporting the pathogenic potential of heterozygous
   variants in patients, albeit with limited penetrance. The relatively
   low phenotypic incidence implies that these variants may act with
   moderate effect sizes and are insufficient to drive disease alone in a
   dominant manner.

   To further investigate the potential additive effects of multiple
   deleterious variants, our analysis of the human cohort revealed an
   elevated TOF risk among individuals with more than two deleterious
   variants in ciliary genes (OR, 1.672; P = 0.0104). This risk increased
   further with variant burden, exceeding a threefold rise when six
   deleterious variants were present, suggesting a cumulative genetic
   effect ([171]Fig. 2). Complementing our human data findings, we then
   generated five double-heterozygous mouse lines by pairwise crossing of
   ciliary gene knockout/mutant strains. Cardiac phenotyping showed a
   33.6% overall incidence of cardiac defects, including OA in 11.48% (14
   of 122) of cases. Notably, the frequency of OA in double-heterozygous
   mice was consistently higher than in single-heterozygous strains, and
   in most cases, also exceeded that observed in homozygous mutants,
   indicating possible synergistic or cumulative effects between different
   ciliary gene variants ([172]Fig. 4). Together, our findings from both
   human genetic analysis and mouse models provide strong evidence for a
   multigenic inheritance mechanism in TOF, in which combinations of
   heterozygous variants across multiple ciliary genes contribute to
   disease risk. The incomplete penetrance observed in heterozygous mice
   further underscores the importance of genetic interactions and
   cumulative variant burden in modulating phenotypic expression.

   Furthermore, beyond the influence of multigenic inheritance pattern,
   the observation of additional cardiac phenotypes such as MU-VSD and ASD
   in ciliary gene knockout/mutant mice may reflect the broader role of
   cilia in cardiac development. Given the widespread distribution of
   cilia across multiple cardiac regions (fig. S4), their dysfunction may
   affect not only OFT-specific processes, resulting in features like OA
   and PM-VSD but also affect other regions involved in septation,
   contributing to MU-VSD and ASD. These phenotypic variations may also be
   influenced by shared molecular pathways among different CHD subtypes
   ([173]13, [174]48), effects of genetic background differences between
   humans and mice ([175]49, [176]50), or environmental factors during
   intrauterine and postnatal development.

   During heart development, cilia serve as transduction centers for
   various development-related signaling pathways, such as the Hh
   signaling pathway, WNT/β-catenin signaling pathway, and TGF-β/BMP
   signaling pathway ([177]51). Previous studies revealed that abnormal
   expression or localization of ciliary genes affects structural
   stability and function of cilia ([178]31). To investigate the common
   mechanisms of cilia in OFT development, we performed RNA-seq of E10.5
   cardiac tissues from wild-type and littermate knockout/mutant mice. The
   DEGs were enriched in pathways related to cell-cell communication
   (including cell adhesion molecules, tight junctions, and ECM-receptor
   interactions), cardiac muscle contraction pathway and Hh signaling
   pathway. Altered activity of cardiomyocyte contraction pathway may
   indicate abnormal cardiac function in knockout and mutant mice. Cell
   adhesion molecules are known to mediate adhesion between cells (mostly
   based on E-cadherin and N-cadherin) and between cells and the ECM
   (integrins). Among them, E-cadherin modulates OFT development by
   activating the Wnt/β-catenin signaling pathway and influencing ISL1
   expression ([179]52, [180]53), whereas integrins regulate OFT
   separation mainly by influencing neural crest cell migration,
   proliferation, and BMP expression ([181]54).

   Mechanistically, we observed shortened cilia length in all six ciliary
   gene knockout and mutant mouse models. Complementary RNA-seq and
   single-cell transcriptomic analyses revealed that ciliary dysfunction
   specifically disrupts Hh pathway activation in cardiomyocytes,
   accompanied by reduced expression of the key transcription factors GLI1
   and GLI2. Previous studies have shown that mutations in cilia-related
   proteins might impair Hh signaling and/or Gli processing, leading to
   severe developmental disorders, while Hh signaling pathway is critical
   for cardiac OFT development and septation ([182]34, [183]35). Notably,
   Gata4 and Nkx2.5 emerged as the only two genes both consistently
   down-regulated across all cilia knockout/mutant strains and identified
   as direct Gli-target genes in the heart by CUT&Tag profiling. As
   expected, the protein expression levels of GATA4 and NKX2-5 were
   prominently decreased in knockout and mutant hearts. GATA4 and Nkx2.5
   act as essential cardiogenic transcriptional regulators implicated in
   many aspects of cardiac development and function. Human genetic studies
   have implicated haploinsufficiency of GATA4 ([184]55, [185]56) and
   Nkx2.5 ([186]57, [187]58) in human TOF. Zhou et al. also reported that
   Gata4 and Gli1 interact genetically, with SHF-specific constitutive Hh
   signaling remarkably rescued CHDs in Gata4 mutant embryos ([188]59).

   Sonic Hh (Shh) signaling has previously been shown to regulate OFT
   alignment by promoting the proliferation of anterior SHF (aSHF)
   progenitors, facilitating myocardial differentiation, and coordinating
   endocardial cushion remodeling ([189]34, [190]40, [191]60, [192]61).
   Shh also maintains Nkx2-5 expression in SHF-derived cardiomyocytes and
   promotes myocardial migration and ECM remodeling within the OFT
   cushions, processes essential for OFT elongation and septation. In line
   with these roles, our data demonstrate that ciliary gene disruption not
   only leads to impaired Hh-Gli signaling and reduced Gata4/Nkx2-5
   expression but also results in shortened OFTs due to reduced
   SHF-derived cardiomyocyte recruitment at E10.5, closely resembling the
   phenotype of Shh/Smo mutant embryos. Our findings extend the classical
   Shh-OFT paradigm by revealing that cilia act as hubs for Hh signal
   amplification in cardiomyocytes, ensuring Gli-dependent transcriptional
   activation of Gata4 and Nkx2-5. This unique mechanism explains how
   subtle perturbations in ciliary function (e.g., mutations) could
   disrupt OFT alignment independent of canonical Shh ligand availability,
   a concept not previously emphasized in developmental studies. While
   prior research has established the essential role of Shh signaling in
   OFT development, our study uncovers an additional regulatory layer in
   which ciliary genes modulate Hh signaling efficiency in cardiomyocytes,
   offering a potential explanation for OFT defects in ciliopathies.
   Moreover, we provide evidence linking Gli1/2 with SHF-derived
   Gata4/Nkx2-5 in heart, thus bridging cilia biology with cardiac
   morphogenesis in a previously unrecognized manner. These insights not
   only enhance our understanding of OFT development but also underscore
   the broader importance of ciliary dysfunction in TOF pathogenesis.

   Although our bioinformatic analyses and mouse models supported the
   pathogenic role of ciliary gene variants in nonsyndromic TOF, this
   study was restricted by several shortcomings. Through our screening
   procedure and enrichment analysis, we found 108 ciliary genes with rare
   LoF variants in our patients with TOF. However, because of the
   polygenic genetic architecture and low penetrance of nonsyndromic TOF,
   it was impossible to construct mouse models for each candidate ciliary
   gene to verify its pathogenic role in TOF. We therefore selected six
   previously unidentified genes (namely AGBL2, MLF1, PKHD1L1, PPEF2,
   SPTBN5, and TEKT3) with recurrent mutations for mouse models
   construction and phenotype validation, both to illustrate the
   reliability of our screening strategy and to provide in vivo evidence
   supporting the contribution of ciliary gene variants to TOF
   pathogenesis. Besides, because of the incomplete penetrance of cardiac
   phenotypes in mouse models, we used homozygous mice for subsequent
   mechanistic studies, as they offer an advantage over heterozygous mice
   in exploring the potential downstream mechanisms of ciliary gene
   mutations affecting TOF pathogenesis.

   In addition, while ciliary dysfunction is known to contribute to a
   spectrum of systemic disorders, this study focused primarily on cardiac
   phenotypes and left-right patterning anomalies (no lateral defects or
   visceral heterotaxy were observed in our models). Whether ciliary genes
   deficiency caused additional abnormal phenotypes in other tissues needs
   to be further studied. Last, this study may still suffer from its
   relatively small sample size due to the challenges in obtaining TOF
   cases, and the absence of parental DNA samples precluded the
   determination of whether the variants were de novo or inherited from
   the unaffected parents. Furthermore, because of ethical and practical
   constraints and given the difficulty of obtaining cardiac tissue from
   living individuals, genetic analyses were based on blood samples rather
   than cardiac tissues. However, as blood may not fully reflect the
   genetic landscape of native cardiomyocytes, future studies
   incorporating cardiac tissue from patients with TOF and matched healthy
   controls will be valuable for uncovering cardiac-specific variants and
   further elucidating the molecular mechanisms underlying TOF.

   In summary, nonsyndromic TOF is a heterogeneous disease with complex
   and diverse etiologies. Our data from next-generation sequencing in 131
   patients systematically identified a higher burden of cilium pathway
   and ciliary gene deleterious variants among individuals with TOF and
   also suggested a multigenic inheritance pattern of ciliary gene
   variants on TOF. Functional studies using mouse models harboring these
   variants recapitulated TOF-like phenotypes as well as demonstrated
   impaired cilia structure and function. Notably, double-heterozygous
   mice exhibited a higher incidence of OA than single-heterozygous
   mutants, further supporting the additive effects of multiple ciliary
   gene variants. Together, our findings highlight the essential role of
   cilia in cardiac OFT development and provide valuable genetic and
   mechanistic insights into the etiology and pathogenesis of nonsyndromic
   TOF. Expanded genetic screening of ciliary genes in larger TOF cohorts
   will be valuable for advancing both diagnosis and precision medicine in
   this population.

MATERIALS AND METHODS

Human sample collection

Patients

   This study enrolled 131 patients who received a confirmed diagnosis of
   nonsyndromic TOF from the First Affiliated Hospital of Nanjing Medical
   University (Nanjing, China) and the Affiliated Children’s Hospital of
   Nanjing Medical University (Nanjing, China) from 2008 to 2014. The
   nonsyndromic TOF samples were diagnosed consistently by pediatric
   cardiologists through comprehensive physical examination,
   echocardiography, and in some cases, further confirmed by surgery.
   Patients were excluded if they showed (i) any extracardiac malformation
   or other types of developmental abnormalities; (ii) chromosomal
   abnormalities; (iii) a positive family history of CHD in their
   first-degree relatives (parent or sibling); (iv) a history of maternal
   diabetes mellitus, phenylketonuria, teratogen exposure, or therapeutic
   drug exposure during pregnancy. For each participant, 5 ml of the whole
   blood were collected for DNA extraction and genomic analyses.

Controls

   For association tests, we applied WGS data of 6004 unrelated
   individuals with noncardiac defects from our Nanjing Lung Cancer Cohort
   project ([193]16) as population genome controls in this study, for
   comparison of gene-based and gene set–based variant frequencies. This
   project was approved by the Human Ethics Committee of Nanjing Medical
   University. Notably, the control group consisted of adults rather than
   children, which was a deliberate choice to mitigate the potential
   underdiagnosis of CHD in younger populations and to reduce selection
   bias in the sample. The sequencing data were processed using the same
   data processing pipeline, and variant annotation were performed in a
   unified manner across case and control samples to minimize potential
   systematic biases due to technical factors or analytical methods.

   Methods and experimental protocols involving human subjects were
   approved by the institutional review boards of Nanjing Medical
   University (2022SCJ2325). All participants completed written informed
   consent form before taking part in this research, and all
   investigations conformed to the principles outlined in the Helsinki
   Declaration.

Animal studies

   All animal experiments were performed under Guidelines for the Care and
   Use of Laboratory Animals and were approved by the Institutional Animal
   Care and Ethical Committee of Nanjing Medical University
   (IACUC-1912049). Mice (C57BL/6J) were housed at a constant temperature
   of 23° ± 2°C, humidity of 55 ± 5%, and under a regular 12-hour
   light-dark cycle in the Animal Research Center of Nanjing Medical
   University. A vaginal plug was checked on the morning following mating
   with an experienced male, and gestational ages were determined by
   visualizing the presence of a vaginal plug (E0.5 = vaginal plug day).
   After mating, pregnant females were individually housed. Before cardiac
   resection, the mice were anesthetized through intraperitoneal
   administration of tribromoethanol (Nanjing Aibei Biotechnology Co. LTD,
   M2920, 0.02 ml/g per mouse), a noncontrolled fast-acting anesthetic
   specifically used for mice, which also allowed for precise dosage
   control. Afterward, the mice were euthanized by cervical dislocation
   and the hearts or embryonic hearts were harvested using a
   stereomicroscope for subsequent experiments.

   All measurements, including PA widths, OFT lengths, cilia lengths and
   proportions, along with histological analysis and CHD diagnosis, were
   conducted by two independent individuals blinded to the mouse genotype.
   In cases of any discrepancies, a third reviewer reevaluated and
   finalized the measurements.

Next-generation sequencing

   WES was performed on 104 patients with nonsyndromic TOF and WGS was
   performed on 27 patients with nonsyndromic TOF. Genomic DNA was
   isolated from blood samples of the subjects using a whole-blood DNA
   purification kit (Qiagen, Hilden, Germany). For WES, exome capture was
   carried out using the Agilent SureSelect Human All Exon V6 Library Prep
   Kit (Agilent, CA, USA), followed by SureSelectXT Target Enrichment
   System for library preparation (Agilent, CA, USA). For WGS, libraries
   were constructed using the standard protocols provided by the
   manufacturer for NEBNext DNA Sample Prep Master Mix Set (NEB, MA, USA).
   All libraries were sequenced with a fragment of 2 × 150–base pair (bp)
   (PE 150) using Illumina HiSeq X Ten sequencers to a depth of 100× (WES)
   and 30× (WGS) base coverage. We used CASAVA (v1.8.2) to convert
   Illumina BCL files to FASTQ files. Then, adapters and low-quality bases
   were trimmed with Trimmomatic (v0.32) ([194]62). After that, sequence
   reads were mapped to the GRCh37 human reference sequence using the
   Burrows-Wheeler Aligner (BWA-MEM v0.7.15-r1140) ([195]63) with the
   default parameters, and the duplicates were marked and discarded using
   SAMBLASTER (v0.1.24) ([196]64). The quality metrics (including
   coverage, median insert size, and percentage of chimeric reads) and
   contamination estimation of “Bam” files were obtained by Picard (v1.70)
   and VerifyBamID2 (v1.04), respectively. Haplotype calling was performed
   using HaplotypeCaller in GATK (v3.8) in GVCF mode according to the best
   practice ([197]65). Individual gVCF files were jointly genotyped for
   high-confidence alleles using Genotype GVCFs for single-nucleotide
   polymorphisms (SNPs) and small insertions and deletions (INDELs) in the
   coding sequences. We performed relatedness and sex checking based on a
   temporary list of high-confidence variants: biallelic, high-call rate
   (>0.99), and common SNPs (allele frequency > 1%). We estimated pairwise
   relatedness by PLINK (v1.9) ([198]66) and inferred the chromosomal sex
   of the samples based on the inbreeding coefficient (F) for common
   variants on chromosome X, excluding pseudo-autosomal regions.

Quality control of sequencing data

Sample quality control

   We excluded individuals with (i) median insert size <250 bp, (ii)
   excess chimeric reads: pct_chimera >0.05, (iii) high contamination:
   freemix >0.05, (iv) unmatched relatedness, and (v) inconsistent
   chromosomal sex and questionnaire information. None of the patients
   with TOF were excluded according to the above criteria, and we finally
   included all of the 131 nonsyndromic TOF cases for further analyses.

Variant quality control

   GATK variant quality score recalibration (VQSR) was used to filter
   variants. The SNP VQSR model was trained using HapMap 3.3 and 1000
   Genomes Omni 2.5 SNP sites, and a 99.6% sensitivity threshold was
   applied to filter variants. We excluded variants if they (i) were on
   sex chromosomes, (ii) failed to pass VQSR filtering criteria, (iii) had
   inbreeding coefficients <−0.3, (vi) had proportions of samples with
   high genotype quality (≥20) ≤0.8, (v) had proportions of samples with
   enough sequencing depth (≥10) ≤0.8, (vi) had call rates ≤0.95, or (vii)
   indels in segment duplication or simple sequence repeat regions.

Variants annotation and identification of potential pathogenic variants in
nonsyndromic TOF

Frequency annotation and the definition of rare variants

   Variant allele frequencies were obtained from the gnomAD (v2.1.1) and
   the ExAC (v0.3.1), both downloaded from
   [199]https://gnomad.broadinstitute.org/downloads, as well as from our
   in-house reference dataset comprising 6004 individuals without CHDs.
   These frequency data were annotated to the variants identified in this
   study using vcfanno (v0.3.2) ([200]67). Variants with estimated MAF_all
   <0.5% in our study cohorts (mutation rate of patients with TOF plus
   control group) and <0.5% in both of the populations [ALL and East Asian
   (EAS)] of gnomAD and ExAC were used for the following analyses.

Definition of potential pathogenic variants of TOF: LoF variants (frameshift,
stop gain, and splice acceptor/donor)

   All coding variants were annotated with SnpEff version 4.3-3 and then
   loaded into the GEMINI software version 0.19.1. For gene set–based
   analyses, all genes annotated to the potential pathogenic variants of
   TOF were used for the following pathway enrichment analysis based on GO
   databases by the ToppGene Suite ([201]https://toppgene.cchmc.org/).

Generation of gene-targeted mouse models

   Ciliary gene knockout and mutant mice were generated by the CRISPR-Cas9
   genome editing system. The single guide RNA (sgRNA) sequence derived
   from each gene (GRCm38/mm10) was designed using CRISPR Design
   ([202]https://crispor.gi.ucsc.edu/), and a DNA oligo corresponding to
   this sequence was cloned into PGL3 for guide RNA in vitro
   transcription. Two single-stranded DNA (ssDNA) homology-directed repair
   (HDR) templates were designed for anthropomorphic variants in Mlf1 and
   Ppef2. A system pool of Cas9 mRNA (25 ng/μl), sgRNA1 (15 ng/μl), and
   sgRNA2 (15 ng/μl) were mixed and microinjected into single-cell mouse
   embryos of the C57BL/6 strain (Agbl2, Pkhd1l1, Sptbn5, and Tekt3); a
   pooled system of Cas9 mRNA (20 ng/μl), sgRNA (20 ng/μl), and ssDNA (1
   μM) conducted the similar experiment as above (Mlf1 and Ppef2). The
   edited sites were amplified by polymerase chain reaction (PCR) and
   Sanger sequencing in the founders, and the founders were followed by at
   least three consecutive backcrosses with congenic mice (C57BL/6) to
   obtain a pure background. All sgRNAs and primer sequences for
   genotyping for each knockout/mutant mouse line are listed in table S4.

Off-target analysis

   The potential OTSs of the sgRNA were predicted with the online CRISPR
   Design Tool ([203]https://cm.jefferson.edu/Off-Spotter/). The top five
   OTSs for each sgRNA were selected, and then the genomic regions
   covering the OTSs were PCR amplified using the primers listed in table
   S6.

Immunofluorescence staining

   To visualize OFT cilia, wild-type and knockout/mutant embryonic hearts
   were harvested at E10.5, E12.5, E14.5, E16.5, and E18.5. Hearts were
   fixed in 4% buffered paraformaldehyde (PFA) (Biosharp, BL539A) on ice
   for 2 hours, rinsed in phosphate-buffered saline (PBS) and incubated in
   20 and 30% sucrose overnight at 4°C until the hearts sunk completely.
   Then, the tissues were embedded in optimal cutting temperature compound
   and stored at −80°C. Cryosections of 15-μm thickness were collected on
   positively charged slides and stored at −20°C. For immunofluorescence
   staining, the hearts were washed in PBS and blocked with donkey serum
   for 1 hour at room temperature before being incubated overnight at 4°C
   with primary antibodies. After three washes with PBS for 5 min, the
   samples were stained for 1 hour at room temperature with secondary
   antibodies followed by 4′,6-diamidino-2-phenylindole (DAPI) staining
   for nuclei visualization. After a second series of washes, the slides
   were mounted in Fluoromount-G (Southern Biotech, 0100-01). Images were
   acquired by a confocal laser scanning microscope (Zeiss, LSM800). The
   antibodies and dilutions were as follows: acetylated tubulin
   (Sigma-Aldrich, T-6793, RRID:AB_477585, 1:1000 dilution),
   acetyl–α-tubulin (Cell Signaling Technology, D20G3, RRID:AB_10544694,
   1:1000 dilution), γ-tubulin (Abcam, ab11317, RRID:AB_297921, 1:1000
   dilution), γ-tubulin (Sigma-Aldrich, T6557, RRID:AB_477584, 1:1000
   dilution), RSPH9 (Proteintech, 23253-1-AP, RRID:AB_2879240, 1:200
   dilution), cardiac troponin t (Abcam, ab8295, RRID:AB_306445, 1:200
   dilution), donkey anti-mouse immunoglobulin G (IgG) Alexa Fluor 488
   (Thermo Fisher Scientific, A21202, RRID:AB_141607, 1:1000), and donkey
   anti-rabbit IgG Alexa Fluor 546 (Thermo Fisher Scientific, A10040,
   RRID:AB_2534016, 1:1000).

Transmission electron tomography

   For TEM, embryonic mouse hearts (E13.5) were excised after cervical
   dislocation and cardiac OFTs were dissected and isolated under a
   stereomicroscope. Then, OFTs were fixed at 4°C overnight with 2.5%
   glutaraldehyde in 0.2 M cacodylate buffer. After washing in cacodylate
   buffer, the OFTs were postfixed in 1% OsO[4] in 0.2 M cacodylate buffer
   for 2 hours at 4°C. Next, the OFTs were washed and poststained with 2%
   aqueous uranyl acetate for 2 hours, dehydrated in a graded ethanol
   series, immersed in acetone and embedded in resin (EPON, Low Viscosity
   Embedding Media Spurr’s Kit #14300). Ultrathin sections (60 to 80 nm)
   were cut on an ultramicrotome (Leica, EM UC7) and mounted on copper
   grids. The sections were then stained in a 2% uranyl acetate saturated
   alcohol solution and lead citrate for 10 min. Images were acquired
   using an FEI transmission electron microscope (FEI, Tecnai G2 Spirit
   Bio TWIN).

Protein extraction and Western blot analysis

   Embryonic mouse heart (E10.5) tissues were collected and lysed using
   universal protein extraction lysis buffer radioimmunoprecipitation
   assay (Beyotime, P0013B) containing a protease inhibitor cocktail
   (MedChemExpress, HY-K0010). The protein lysates were centrifuged and
   the supernatants were quantified using the BCA protein assay kit
   (Beyotime, P0012). Then, the lysates were denatured at 100°C for 10 min
   with 5× SDS–polyacrylamide gel electrophoresis sample buffer (Beyotime,
   P0015L). For Western blotting, 40 μg of protein was separated on 8% or
   10% SurePAGE according to the molecular weight of the protein
   (GenScript, M00661/M00666) and transferred onto a polyvinylidene
   difluoride membrane (Millipore, ISEQ00010). The antibodies used in
   Western blotting are as follows: AGBL2 (Invitrogen, PA5-101547,
   RRID:AB_2850981, 1:1000 dilution), MLF1 (Santa Cruz Biotechnology,
   sc-514294, 1:500 dilution), PPEF2 (Abcam, ab183667, 1:1000 dilution),
   TEKT3 (Invitrogen, PA5-70344, RRID:AB_2690235, 1:1000 dilution), GLI1
   (Proteintech, 66905-1-Ig, RRID:AB_2882232, 1:1000 dilution), GLI2 (R&D
   Systems, AF3635, RRID:AB_2111902, 1:1000 dilution), GATA4 (Santa Cruz
   Biotechnology, sc-25310, RRID:AB_627667, 1:500 dilution), NKX2-5
   (Abcam, ab97355, RRID:AB_10680260, 1:1000 dilution), and GAPDH
   (glyceraldehyde-3-phosphate dehydrogenase; Cell Signaling Technology,
   5174, RRID:AB_10622025, 1:2000 dilution). Knockout/mutation of AGBL2,
   MLF1, PPEF2 and TEKT3 were verified by Western blotting, whereas
   validation of PKHD1L1 or SPTBN5 knockout in mice was performed using
   immunofluorescence with a PKHD1L1 (Abnova, K7241-1F5, 1:1000 dilution)
   or SPTBN5 antibody (Abcam, ab111223, RRID:AB_10864996, 1:1000 dilution)
   respectively, as described above.

Histological analysis of mouse hearts

   For histological examination, the hearts of wild-type and
   knockout/mutant newborns at P1 were harvested and fixed overnight in 4%
   PFA (Biosharp, BL539A) and then dehydrated in a graded ethanol series,
   cleared, and paraffin wax embedded. Transverse sections (6 μm) were
   stained with H&E using standard histology techniques. The stained
   hearts were photographed with a stereomicroscope (Nikon, SMZ1270).

MEF isolation and location analysis

   Primary MEFs were isolated from E13.5 embryos. Briefly, embryos were
   minced using a sterile razor blade. Then, 4 ml of 0.25% trypsin-EDTA
   (Gibco, 25200-056) was added, and the samples were incubated at 37°C
   for 15 min. After incubation for 15 min in 10-cm tissue culture dishes,
   the remaining fragments were pipetted multiple times to dissociate the
   cells and incubated overnight in Dulbecco’s modified Eagle’s
   medium/Nutrient mixture F-12 (DMEM/F-12; Gibco, C11330500BT)
   supplemented with 10% fetal calf serum (Gibco, 10099141C). For primary
   cilia analysis, MEFs were serum starved in DMEM/F-12 supplemented with
   0.1% serum for 24 hours to induce cilia growth. Cells were fixed with
   4% PFA in PBS and immunostained with antibodies against AGBL2
   (Invitrogen, PA5-101547, RRID:AB_2850981, 1:200 dilution), MLF1 (Santa
   Cruz Biotechnology, sc-514294, 1:200 dilution), PPEF2 (Abcam, ab183667,
   1:200 dilution), PKHD1L1 (Abnova, K7241-1F5, 1:200 dilution), SPTBN5
   (Abcam, ab111223, RRID:AB_10864996, 1:200 dilution), and TEKT3
   (Invitrogen, PA5-70344, RRID:AB_2690235, 1:200 dilution), together with
   structural markers acetylated tubulin (Sigma-Aldrich, T-6793,
   RRID:AB_477585, 1:1000 dilution), acetyl–α-tubulin (Cell Signaling
   Technology, D20G3, RRID:AB_10544694, 1:1000 dilution), γ-tubulin
   (Abcam, ab11317, RRID:AB_297921, 1:1000 dilution), γ-tubulin
   (Sigma-Aldrich, T6557, RRID:AB_477584, 1:1000 dilution), and Nek2
   (Santa Cruz Biotechnology, sc-55601, RRID:AB_1126558, 1:200 dilution).
   The cells were then imaged using a confocal laser scanning microscope
   (Zeiss, LSM800).

Quantification of cilia

   The percentage of ciliated cells and the length of cilia were
   quantified and measured on 15-μm-thick heart sections using ZEN
   software. Z-stacks were set by finding the highest and lowest depths
   (with 63× objectives) with visible fluorescence and using the system
   optimized setting to determine steps. Z-stacks were then compiled to
   form maximum projection images. Three-dimensional reconstructions of
   these images were performed by importing Z-stack confocal images into
   Imaris 9.0 (Bitplane Inc.) and creating surface renderings based on
   stain intensities. Cilia length was measured from the base of the
   axoneme (acetylated tubulin–positive staining) to the tip. The
   percentage of ciliated cells was calculated as the number of cilia
   (acetylated tubulin and γ-tubulin–positive staining)/number of
   DAPI-labeled nuclei for each field.

RNA-seq and data analysis

   Total RNA from the heart tissues of ciliary gene knockout/mutant mice
   and wild-type littermate mice at E10.5 was extracted with TRIzol X-100
   (Thermo Fisher Scientific, 15596018) and a total amount of 2 μg of RNA
   per sample was used for RNA library construction. RNA integrity and
   concentrations were assessed using the RNA Nano 6000 Assay Kit of the
   Bioanalyzer 2100 system (Agilent Technologies). Sequencing libraries
   were then generated using the NEBNext Ultra RNA Library Prep Kit
   (Illumina, E7530). mRNA was purified from total RNA using poly-T
   oligo–attached magnetic beads and cut into small fragments for cDNA
   synthesis. After adenylation of the 3′ ends of DNA fragments, adaptors
   with hairpin loop structures were ligated, and then PCR was performed.
   High-throughput RNA sequencing was carried out using Illumina
   NovaSeq6000 platform and 150-bp paired-end reads were generated.

   Raw data were filtered and processed to obtain clean data with high
   quality and aligned to the mouse reference genome GRCm38/mm10 with STAR
   software (v2.5.3a). Then, assembly and quantification of the
   transcripts were accomplished with RESM software (v1.3.3) guided by the
   Ensembl gtf gene annotation file
   ([204]https://ftp.ensembl.org/pub/release-93/gtf/mus_musculus/). Read
   counts were used for the measurements of the relative abundance of the
   transcripts. Fold change and P values were calculated on the basis of
   read counts using edgeR package (v3.36.0), and DEGs were identified
   (FDR < 0.05). KEGG pathway enrichment analysis was performed on the
   basis of the DEGs by the ToppGene Suite
   ([205]https://toppgene.cchmc.org/).

scRNA-seq and data analysis

   E10.5 wild-type and homozygous ciliary gene knockout/mutant mouse
   embryos were isolated and dissected for cardiac OFTs. The obtained
   tissues were washed for two to three times with Dulbecco’s
   phosphate-buffered saline (DPBS) precooled at 4°C and then was
   transferred to a 2-ml centrifuge tube with tissue storage buffer
   (Miltenyi Biotec, #130-100-008). The tube was centrifuged at 50g for 1
   min at 4°C, after which the tissue storage buffer was discarded, and 1
   ml of prewarmed digestion solution (Miltenyi Biotec, #130-110-203;
   37°C) was added to the OFT tissue. The tissues were incubated in a 37°C
   water bath and gently pipetted every 5 min, for a total digestion time
   of approximately 15 min. Enzymatic digestion was terminated by adding
   10% fetal bovine serum (Pricella, #164210-50), and the cell suspension
   was filtered through a 40-μm cell strainer (Corning, #CLS431750). The
   filtered suspension was centrifuged at 300g for 3 min at 4°C, and the
   pellets were resuspended in PBS containing 0.01% bovine serum albumin
   (Sigma-Aldrich, #A1933). After treating with AO/PI double fluorescent
   dye, cell quality was conducted on Countstar Rigel S2 to determine cell
   concentration, cell viability, and cell clump. At the same time, Trypan
   blue staining was performed to observe the cell morphology, cell
   concentration, cell fragments and impurities.

   Single-cell encapsulation was performed using the Chip A Single Cell
   Kit v2.0 (MobiDrop, #S050100201) on the MobiNova-100 platform
   (MobiDrop, #A1A40001). Each droplet contained one cell and a gel bead
   carrying millions of oligonucleotides with unique cell barcodes.
   Droplets were then processed with MobiNovaSP-100 (MobiDrop, #A2A40001),
   which initiated oligo release and hybridization with cellular mRNAs.
   Reverse transcription and cDNA amplification were performed within
   droplets. Libraries were constructed using the High-Throughput Single
   Cell 3′RNA-Seq Kit v2.0 (MobiDrop, #S050200201) and indexed using the
   3′ Single Index Kit (MobiDrop, #S050300201). Sequencing was carried out
   on an Illumina NovaSeq 6000 platform using 150-bp paired-end reads.

   Raw sequencing data (FASTQ files) were demultiplexed and aligned to the
   mouse reference genome (GRCm38/mm10) using Cell Ranger (v6.1.1)
   pipeline with default parameters, generating a feature-barcode
   expression matrix. Downstream analysis was performed in R using the
   Seurat package (v4.1.3). Genes expressed in fewer than three cells, as
   well as cells with <200 genes or >10% mitochondrial gene content, were
   excluded. Doublets were identified and removed using the DoubletFinder
   package (v2.0.3). Batch correction was conducted using the Harmony
   package (v0.1.1) to minimize potential sample-specific effects. Eleven
   significant principal components were selected on the basis of
   ElbowPlot and JackStrawPlot outputs in Seurat. Cell clustering and
   dimensionality reduction were performed using the Seurat FindClusters
   and RunUMAP functions, respectively. Cluster-specific marker genes were
   identified using the FindAllMarkers function. Cell types were annotated
   on the basis of known marker gene expression.

Gene set variation analysis

   GSVA was performed to evaluate Hh signaling pathway activity in
   cardiomyocytes from wild-type and ciliary gene knockout/mutant mice
   ([206]68). A curated mouse gene set for Hh signaling
   ([207]HALLMARK_HEDGEHOG_SIGNALING.v2023.1.Mm) was downloaded from the
   GSEA database. This gene set includes representative components of the
   Hh pathway such as GLI family genes, Shh, Gli1, Ptch1, and others
   directly related to Hh signaling. GSVA was conducted using the “GSVA” R
   package (version 1.40.0) with default parameters. Enrichment scores
   were calculated for each cardiomyocyte to represent the relative
   activity level of the Hh signaling pathway based on gene expression.

Cleavage under targets and tagmentation

   CUT&Tag assays and subsequent high-throughput sequencing were performed
   by Seqhealth Technology Co. Ltd. (Wuhan, China). For each group, 10
   fresh E10.5 hearts from wild-type mouse embryos were collected and
   homogenized to isolate cells, which were then incubated with
   concanavalin A–coated magnetic beads. The cells were subsequently
   incubated with primary antibodies against GLI1 and GLI2 (GLI1: R&D
   Systems, AF3324, RRID:AB_663902; GLI2: R&D Systems, AF3526,
   RRID:AB_2279108) for 2 hours, followed by incubation with a secondary
   antibody for 1 hour at room temperature. Hyperactive pG-Tn5 transposase
   was added to initiate tagmentation. After terminating the reaction, DNA
   fragments were purified using PCI extraction and amplified via PCR with
   indexed P5 and P7 primers. Final libraries were enriched, quantified,
   and sequenced on an Illumina NovaSeq 6000 platform using a PE150
   configuration.

   Raw sequencing data were processed using Fastp (v0.23.1) for quality
   filtering and adapter trimming. Clean reads were aligned to the
   GRCm38/mm10 mouse reference genome using Bowtie2 (v2.2.6) with default
   settings. Sambamba (v0.7.1) was used to convert file formats and remove
   PCR duplicates. DeepTools (v2.4.1) was used to visualize read
   distributions around transcription start sites. Peak calling was
   conducted using MACS2 (v2.2.7.1), and peak annotation and distribution
   analysis were performed using bedtools (v2.30.0). Motif analysis was
   carried out using HOMER (v4.10). Last, transcription factor binding
   sites were predicted on the basis of enriched peak sequences and motif
   data from the JASPAR database.

Statistical analysis

   Gene-based association analyses were performed by aggregating the
   qualifying rare deleterious variants into genes and comparing allele
   counts between TOF cases and controls using two-sided Fisher’s exact
   tests. Benjamini-Hochberg correction was applied, with an FDR threshold
   of <0.2. Input data consisted of rare, predicted deleterious variants
   (e.g., LoF and detrimental missense variants) and filtered on the basis
   of MAF thresholds (<0.005) and functional annotations. The LoF model
   focused on LoF variants, defined as frameshift, stop-gain, and
   canonical splice site (splice acceptor/donor) variants. The detrimental
   missense variant model applied four prediction criteria: (i) predicted
   as pathogenic by AlphaMissense (score, > 0.564) ([208]69); (ii) CADD
   score, >20 ([209]70); (iii) predicted as deleterious by SIFT
   (score, ≤ 0.05) ([210]71); and (iv) predicted as probably_damaging
   (score, 0.447 to 0.909) or possibly_damaging (score, 0.909 to 1) by
   PolyPhen2 ([211]72). Each criterion was applied independently to
   identify genes enriched for detrimental missense variants in TOF cases.
   Gene set–based association analyses were conducted for LoF variants and
   detrimental missense variants in all ciliary genes involved in cilium
   pathway between TOF cases and controls using two-sided Fisher’s exact
   test. Missense variants were considered deleterious only if all four
   criteria were simultaneously met to ensure high-confidence
   classification in oligogenic burden assessment. ORs were calculated to
   evaluate the burden of varying numbers of deleterious variants in
   ciliary genes among patients with TOF compared to control individuals,
   with statistical significance defined as P < 0.05 and OR > 1.

   For multigenic burden analysis, variant counts per individual were
   based on the number of unique affected genes. Multiple variants within
   the same gene were consolidated and counted once, whereas variants in
   different genes were counted individually on the basis of their actual
   occurrence. This approach ensured that the per-sample mutation burden
   reflected the number of distinct genes affected.

   For in vitro and in vivo experiments, the results were presented as the
   means ± SEM. Statistical calculations were performed using unpaired
   two-tailed Student’s t test or Wilcoxon–Mann-Whitney test for
   continuous variables data and Pearson’s Chi-squared test or Fisher’s
   exact test for categorical data, as appropriate. A value of P < 0.05
   (*) was considered statistically significant, with **P < 0.01 and
   ***P < 0.001.

Acknowledgments