The transcriptomic datasets were retrieved from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/, including GSE203274 for single-nucleic RNA-seq, GSE146220 for a microarray dataset, and GSE217772 for bulk RNA-seq). GSE203274 includes the Fastq form of single-nucleic RNA-seq data based on Illumina Novaseq 6000 platform (15). Six cardiac tissues of TOF (SRR19266878, SRR19266879, SRR19266880, SRR19266881, SRR19266882, SRR19266883) and three tissue of normal (SRR19266894, SRR19266896, SRR19266898) were utilized in our study. Two bulk expression profiles were also utilized. GSE146220 is a microarray dataset containing five heart tissues of patients diagnosed as Tetralogy of Fallot and four from control. The study has not yet been published. GSE217772 is a bulk RNA-seq dataset based on Illumina NovaSeq 6000. This dataset includes five cardiac tissues of patient with Tetralogy of Fallot and five from heathy control, with the study not being published neither.

The raw data was retrieved and transformed into fastq files using SRA toolkit. These datasets were further transformed to count matrixes using the Cellranger sofware (v7.0.1) based on homo sapiens reference genome GRCh38. The Seurat framework was subsequently employed for the downstream analysis (16). Briefly, a minimum of 500 genes and 35% mitochondrial was set as the cutoff for the removement of low-quality cells because cardiomyocytes contain more mitochondria. The batch effects were eliminated using the Harmony algorithm (17). Harmony achieves data iteration by clustering similar cells from different batches while maximizing the heterogeneity of batches within each cluster. In each iteration, Harmony computes a correction factor for individual cells and applies it to the data, effectively mitigating batch effects while preserving inherent biological variations. Doublets were eliminated using the R package “DoubletFinder” (17). Data preprocessing was based on cell cycle genes and known technical noise features. Then, DoubletFinder simulated dual cell populations to generate artificial double cells, which were then compared with the actual data to identify and eliminate the double cell effect. This approach significantly enhanced the accuracy and reliability of our dataset. The scaled data with top 3,000 variable genes were used to perform principal component analysis (PCA). The first 15 principal components were used for further clustering. To reduce the dimensionality of snRNA-seq data, we employed the Uniform Manifold Approximation and Projection (UMAP) technique. UMAP is a renowned method for dimensionality reduction that effectively preserves both global and local structures of the data in a lower-dimensional space. Marker genes of each type of cell were curated in previous study (15). Differentially expressed genes (DEGs) for each cell cluster between TOF and normal samples were identified using the “FindMarkers” function (Wilcoxon rank-sum test) within the Seurat package, following aggregation into “pseudobulk” samples to enhance statistical power. Genes were considered differentially expressed if they exhibited an |log2FC| > 0.25 and an FDR-adjusted p-value <0.1. Pseudobulk analysis was conducted by aggregating 50 cells within a cluster to form a pseudobulk sample. The total count of a gene was used as gene expression. Expression density of specific genes was exhibited using the R package “Nebulosa” (18). Single sample gene set enrichment analysis (ssGSEA) was conducted to score the expression of genes at single-cell resolution.

The microarray expression matrix was log transformed for further analysis. DEGs were identified from the microarray dataset (GSE146220) using the “limma” R package (19) and from the RNA-seq dataset (GSE217772) using the “edgeR” R package (20). For both analyses, genes with an absolute log2 fold change (|log2FC|) >1 and a false discovery rate (FDR) adjusted p-value <0.05 were considered statistically significant. The complete lists of DEGs from both datasets have been included as Supplementary Tables S1, S2. Functional enrichment analysis was performed using the webtool Metascape (https://metascape.org/gp/index.html) (21). The protein-protein interaction (PPI) network was constructed using the STRING web tool (v11.0). Only interactions with a combined score >0.7 (high confidence) were included in the network. Disconnected nodes were hidden from the final network visualization.

The samples used for qRT-PCR analysis were obtained from the outflow tract portion of fetal TOF heart tissues, collected from women who voluntarily terminated their pregnancies during the second trimester, a critical period for diagnosing congenital heart defects by ultrasound. TOF fetuses with other developmental abnormalities were excluded. Normal control outflow tract samples were obtained from second-trimester therapeutic terminations for severe maternal health reasons. Echocardiography confirmed that these fetuses had normal cardiac development and no other developmental defects. The Medical Ethics Committee of Sichuan Provincial People's Hospital approved the study protocol and the collection of human samples, and written informed consent was obtained from the parents of each participant (No. 016292). All experiments complied with the Helsinki Declaration and relevant national laws.

Total RNA was extracted from heart tissues using TRIzol reagent (Life Technologies, Carlsbad, CA, USA). Reverse transcription of the RNA was performed using a high-capacity cDNA reverse transcription kit (OSAKA, JAPAN). Quantitative real-time PCR (qRT-PCR) was carried out with SYBR Green PCR Master Mix (Roche, Germany) according to the manufacturer's instructions, on a BIO-GENER Real-Time System (China). Each sample was run in triplicate. The cycle threshold (CT) values were normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH) measured on the same plate, and the 2^−ΔΔCT method was used to calculate fold differences in gene expression (22). We selected genes with high node degree in the PPI network for validation in the human TOF group and normal control group. The upregulated genes (UpGs) selected include UQCRC1, UQCRQ, HADHA, and HADHB. The downregulated genes (DpGs) selected include COL5A1, COL1A2, THBS2, TAGLN, and ELN. The primer sequences for the target genes and reference genes have been provided in a table format within the supplementary materials (Supplementary Table S3).

All the statistics was performed using the R software (v4.3.1). Wilcoxon test was employed for the calculation of DEGs of each cluster of snRNA-seq data. The p-value was adjusted by False Discovery Rate (FDR). FDR < 0.1 was considered statistically significant in DEGs identification and enrichment analysis unless otherwise specified.

The bulk data sets of cardiac tissue of TOF and healthy control were first leveraged to explore genes potentially participated in disease progression. DEGs were calculated and exhibited (Figure 1A). Describing the association between genes labeled out in Figure 1B with TOF-related signaling pathways (Figure 1B). The consistent DEGs in two cohorts were intersected, resulting in 29 consistent UpGs (Figure 2B) and 22 DpGs (Figure 2B). were constructed, indicating intense interactions between these genes. We used the STRING network tool to exclude genes that did not have significant interactions with other genes and constructed the PPI networks. In the constructed PPI network, UpGs such as UQCRQ, UQCRC1, HADHA, and HADHB are identified as hub genes, meaning they have the most interaction nodes (Figure 2C). These genes are involved in components of mitochondrial respiratory chain complex III (UQCRQ, UQCRC1) and key enzymes in the fatty acid β-oxidation pathway (HADHA, HADHB). In the constructed PPI network of DpGs, genes such as COL1A2 and COL5A1 are identified as hub genes, meaning they have the highest number of interaction nodes (Figure 2D). These genes are involved in collagen biosynthesis and extracellular matrix (ECM) remodeling. Functional enrichment analysis was performed to delineate the functions of consistent up- or down-regulated genes based four distinct databases. As a result, these DEGs were top enriched in biological process including aerobic electron transport chain, heart morphogenesis, and extracellular matrix organization, indicating an implication in TOF (Figure 3A). The KEGG enrichment analysis of DEGs highlighted significant involvement in pathways such as diabetic cardiomyopathy, fatty acid degradation, protein digestion and absorption, and focal adhesion (Figure 3B). The Reactome enrichment analysis of DEGs revealed significant involvement in pathways such as mitochondrial fatty acid beta-oxidation of unsaturated fatty acids, the citrate (TCA) cycle and respiratory electron transport, mitochondrial translation initiation, and others (Figure 3C). The Wiki enrichment analysis of DEGs highlighted significant involvement in pathways such as fatty acid β-oxidation, the electron transport chain (OXPHOS system in mitochondria), etc. (Figure 3D). Interestingly, fatty-acid beta-oxidation and related pathways were consistently enriched, indicating an association between fatty acid metabolism and TOF.

To further explore the molecular underpinnings of TOF at single-cell resolution, expression profiles of six cardiac tissues of TOF and three normal sample underwent snRNA-seq were collected. Stringent quality control procedures were conducted, including Harmony-based batch effect elimination, DoubletFinder-based doublet detection, providing a reliable foundation for subsequent analysis (Supplementary Figure S1). As a result, a total of 135,646 quantified cells were identified, with 50,911 from normal and 84,735 from TOF (Figures 4A–C). Unbiased clustering identified 11 types of cells, including cardiomyocyte, endothelial cell (endothelial), smooth muscle cells (SMC), cardiofibroblast, epithelial-like cell, pericyte, monocyte/macrophage, lymphocyte, and a small amount of neuron (Figure 4A). The proportion of different cell subpopulations in the cardiac tissue of the two sample groups differs. Compared to the normal control group, the TOF group shows an increased proportion of monocytes/macrophages, epithelial-like cells, cardiomyocytes, smooth muscle cells (SMCs), endothelial cells, and others in the cardiac tissue (Figure 4B). The heatmap displaying differentially expressed genes of each cell type highlights marker genes used to distinguish between various cell subpopulations (Figure 4C). We have defined up-regulated and DpGs that may be involved in the pathology of TOF. UMAP embeddings disclosed that the UpGs were enriched in cardiomyocytes, and the DpGs were enriched in cardiofibroblast (Figure 5A). The TOF tissue exhibited significantly elevated ssGSEA scores in fatty acid metabolism-associated pathways compared to normal controls. Specifically, these pathways included fatty acid metabolism, beta oxidation of very long chain fatty acids, and omega-3 and omega-6 fatty acids synthesis (Figure 5B). These findings collectively indicate that TOF tissue possesses a markedly higher level of fatty acid metabolic activity. Notably, our pathway enrichment analysis revealed that these fatty acid metabolic pathways were predominantly active within cardiomyocytes (Figure 5C). This cell-specific enrichment suggests that dysregulated fatty acid metabolism in cardiomyocytes may contribute significantly to the pathological mechanisms underlying TOF, potentially driving disease progression through impaired energy homeostasis or lipotoxic effects.

To confirm the findings from the bioinformatics analyses, qRT-PCR was used to assess the expression of the hub genes in normal and TOF heart tissues. In line with the bioinformatics predictions, TOF heart tissues showed an increased expression of UQCRC1 (Figure 6A), UQCRQ (Figure 6B), and HADHA (Figure 6C), along with a decreased expression of COL5A1 (Figure 6E), COL1A2 (Figure 6F), TAGLN (Figure 6H), and ELN (Figure 6I) compared to normal tissues. The expression of HADHB (Figure 6D) and THBS2 (Figure 6G) showed no difference between the two groups.