This study received approval from the Institutional Review Board at the Children’s Hospital of Philadelphia (CHOP). PBMCs from three de-identified children (2 males and 1 female) with SV/HLHS were compared to those from three healthy controls (2 males and 1 female). The study design included three biologically independent SV/HLHS cases and three matched healthy controls, providing three replicates per group. This level of replication is commonly accepted in exploratory studies, where the primary goal is to identify robust and biologically meaningful patterns. Blood samples were collected in EDTA-coated tubes and promptly processed at the Center for Applied Genomics (CAG) at CHOP. PBMCs were isolated using Ficoll density gradient centrifugation. scRNA-seq for each sample was conducted using the 10X Chromium Single Cell Gene Expression assay (10x Genomics, Single Cell 3' v3) [10]. Sequencing was performed on the Illumina HiSeq2500 SBS v4 platform. The resulting data from the Chromium scRNA-seq were processed with the Cell Ranger 7.1.0 software suite (10x Genomics), with sequencing reads aligned to the GRCh38 reference genome. The number of cells sequenced and analyzed for each sample ranged from 4,598 to 14,298.

The scRNA-seq data were analyzed using the Seurat R package (v 5.1.0) [11, 12]. Cells were retained if they had more than 200 and fewer than 5,000 detected genes and <15% mitochondrial‐gene reads. Gene expression was normalized with the LogNormalize method (scale factor = 10,000). We identified the top 2,000 variable genes with the variance-stabilizing transformation (vst) method, scaled all genes, and performed principal-component analysis (PCA). The first 10 PCs were used to construct a shared-nearest-neighbor (SNN) graph, and clusters were called with the Louvain algorithm (resolution = 0.5). Uniform manifold approximation and projection (UMAP) based on the same PCs provided two-dimensional visualization. Marker genes were defined using min.pct = 0.25 and |log₂FC| > 0.25.

Natural-killer (NK) cells were identified with SingleR in combination with celldex::DatabaseImmuneCellExpressionData () [13], yielding 548 – 939 NK cells per sample after filtering. To resolve NK-cell heterogeneity, we extracted the normalized NK-cell expression matrices, transposed them (rows = cells, columns = genes), and performed unsupervised K-means clustering with scikit-learn (v 1.4.0) [14]. Candidate cluster numbers k = 2 – 10 were evaluated, and the solution with the highest average silhouette score was chosen; clustering was run with random_state = 0 for full reproducibility. Outliers were defined as cells lying more than 2 SD beyond the mean Euclidean distance to their cluster centroid. Cluster structure was visualized with t-distributed stochastic neighbor embedding (t-SNE) [15], and plots were generated in matplotlib [16]. Cluster assignments and outlier flags are provided as Supplementary Data.

Gene-set enrichment analysis (GSEA v 4.3.2) [17] was carried out against the Hallmark collection [18] of the Molecular Signatures Database (MSigDB) [19]. Differential-expression testing between NK-cell clusters (excluding genes used for K-means clustering) was performed on log-transformed counts with a two-sample, two-tailed Student’s t-test. We used the CellChat R package (version 1.6.1) [17] to infer intercellular communication networks from our normalized NK cell expression data.

In our study, we observed that NK cell clusters tend to be heterogeneous (Figure 1). With their function still unknown, there are no definitive gene markers available to classify these NK subpopulations with conventional scRNA-seq pipelines or traditional flow cytometry. To test whether these heterogeneous NK cells are responsible for the DE genes, we examined the 1600 DE genes identified in NK cells from our previous study [9]. K-means clustering with these DE genes yielded high Silhouette scores with 2 clusters in all the 6 samples (Figure 2). As shown in Figure 3, the NK cells within the same sample can be divided into two distinct clusters. Despite clinical heterogeneity among the cases, including differences in anatomical diagnoses, surgical histories, and comorbid conditions, we observed consistent NK cell clustering patterns and transcriptomic signatures across all three samples, suggesting a potentially conserved NK cell transcriptomic signature in SV/HLHS. We then performed GSEA using the averaged difference between the two clusters across the three cases (Table 1; Supplementary Table S1). In one cluster, the gene sets HALLMARK_ANDROGEN_RESPONSE and HALLMARK_HYPOXIA are significantly upregulated, and HALLMARK_HEME_METABOLISM is significantly downregulated.

Compared to the controls, the upregulated androgen response and downregulated heme metabolism are more significant in the cases, as shown by Cluster expression in Cases vs. Controls in Table 1. In addition, upregulation of the gene sets HALLMARK_OXIDATIVE_PHOSPHORYLATION, HALLMARK_INTERFERON_ALPHA_RESPONSE, and HALLMARK_INTERFERON_GAMMA_RESPONSE are more significant in Cluster 1.

To further characterize the two NK cell clusters, we examined the expression of the two NK cell clusters for genes not used for the clustering. We identified 789 genes with DE between the two clusters in all three cases (Supplementary Table S2). GSEA analysis showed upregulated HALLMARK_ALLOGRAFT_REJECTION and downregulated HALLMARK_UV_RESPONSE_DN in Cluster 1 in cases. The downregulated HALLMARK_UV_RESPONSE_DN in Cluster 1 in cases is more significant than that in controls (Table 1). CellChat-inferred ligand–receptor interactions between Cluster 0 and Cluster 1 are listed in Supplementary Table S3.

In our analysis of NK cell populations, Cluster 1 exhibits distinct transcriptional profiles compared to Cluster 0 in term of the DE genes associated with SV/HLHS, marked by significant upregulation of gene sets including HALLMARK_ANDROGEN_RESPONSE and HALLMARK_HYPOXIA, alongside a notable downregulation of HALLMARK_HEME_METABOLISM. Androgen signaling may modulate NK cell immune functions, potentially tempering their cytotoxic activity while promoting survival and adaptation, particularly under stress [20]. Concurrently, the increased expression of hypoxia-responsive genes likely enhances NK cells' capability to operate in oxygen-deprived environments [21]. The reduced heme metabolism complements the other observed changes by potentially reducing oxidative stress and regulating metabolic responses, including oxidative phosphorylation, amino acid, and xenobiotic metabolism [22], likely shifting towards more energy-conserving processes and potentially linking Cluster 1 to the elevated thrombotic risk observed clinically. Collectively, these transcriptional changes suggest that NK cells in Cluster 1 are likely geared towards adaptation and survival in hypoxic or stress-related conditions. In contrast, Cluster 0 with different gene expression patterns might retain a higher cytotoxic capability and responsiveness, better suited for environments requiring rapid and robust immune reactions without the adaptive pressures of chronic stress.

The distinct overexpression of cytotoxic effector genes NKG7, GNLY, GZMB, and FGFBP2 in Cluster 1 versus key transcription factors ZEB1, BACH2, and BCL11B in the other reflects divergent, subtype-specific regulatory networks and stable epigenetic landscapes. Studies show BACH2 loss drives a discrete cytotoxic program with increased granzymes rather than a gradual shift [23], ZEB1 expression is associated with maintenance/less-cytotoxic states by repressing effector loci [24], and BCL11B drives a distinct adaptive-like differentiation incompatible with the cytotoxic program [25]. In addition, epigenetic analyses reveal fixed chromatin states for each subtype [26], and single-cell clustering consistently uncovers separate NK endpoints with no intermediate transcriptional profiles [27]. These marker genes show differential expression in opposite directions between clusters, emphasizing their divergent functional programs without evidence of a gradual transition pattern. Together, these patterns support the interpretation that the subtypes represent two functionally opposite, terminal NK cell states rather than stages along a developmental or activation trajectory.

In addition, our cell–cell communication analysis revealed that the cytotoxic cluster (Cluster 1) preferentially engages inhibitory checkpoints such as HLA-E binding to the CD94/NKG2A heterodimer [28, 29] and TGFB1 signaling through TGFβR1/2 [30], suggesting built-in mechanisms to restrain potent effector activity, whereas the transcription factor–high cluster (Cluster 0) shows alternative contact-mediated interactions (e.g., CLEC2–KLRB1) [31], consistent with niche or maintenance signals. These divergent communication patterns align with the distinct overexpression of the DE genes described above. Together, they reinforce that each subtype represents a terminal NK cell state with its own regulatory network and signaling milieu.

The two clusters of NK cells observed in our study, present in both healthy controls and cases, likely represent distinct functional states tailored to specific immunological needs. Notably, marked differences in the expression profiles of these clusters between cases and controls were observed, as shown by Cluster expression in Cases vs. Controls in Table 1. Cluster 1 in the cases shows significantly heightened androgen response and more pronounced downregulation of heme metabolism compared to controls. These changes may serve a dual role. This modulation could be protective, reducing potential tissue damage in chronic disease contexts such as seen in SV/HLHS. On the other hand, this reduced cytotoxicity might compromise the NK cells' ability to clear pathogens effectively, potentially increasing the risk of infections. In addition, upregulation of the gene sets related to oxidative phosphorylation, and interferon responses, both alpha and gamma, are more significant in Cluster 1 in cases. Enhanced oxidative phosphorylation indicates that these NK cells have elevated energy production capabilities [32]. Upregulated Interferon responses imply this cluster of NK cells is primed for responding to viral infections and potentially other pathogens [33]. Altogether, upregulated HALLMARK_ANDROGEN_RESPONSE and downregulated HALLMARK_HEME_METABOLISM, plus these additional adaptations, suggest strategic modifications of NK cells to meet the demands of specific pathological states in prolonged exposure to pathogens or inflammatory conditions in SV/HLHS. These findings raise the possibility that pharmacologic modulation of these pathways could help restore a more balanced NK-cell phenotype, enhancing cytotoxicity while preserving protective adaptation in SV/HLHS.

In our extended analysis, we examined additional gene expressions within the two NK cell clusters beyond those initially used for clustering. This examination revealed 789 genes with DE between the two clusters across all cases. Specifically, Cluster 1 in cases exhibited upregulated HALLMARK_ALLOGRAFT_REJECTION and downregulated HALLMARK_UV_RESPONSE_DN. The upregulation of the allograft rejection gene set in Cluster 1 underscores a heightened immunological readiness to recognize and attack non-self cells, a crucial feature for defending against foreign tissues and potentially harmful pathogens [34]. The downregulation of genes associated with the UV response, including key regulators like DYRK1A, ATXN1, ATP2B1, NIPBL, RUNX1, and SIPA1L1, presents a complex scenario. These genes are crucial for various NK cell functions such as signaling, development, and effector responses [35, 36]. Their reduced expression might impair the NK cells’ capacity for regular immune surveillance and cellular migration. Furthermore, the downregulation of genes associated with the UV response in cases is significantly higher than controls. When considered together, these opposing trends in gene expression paint a picture of a dual nature in NK cells in Cluster 1, which could render these NK cells less effective in handling routine immune tasks, posing risks for viral infections and complicate the management of SV/HLHS.

A focused transcriptomic signature derived from these findings could be developed for longitudinal monitoring, with potential to correlate expression dynamics with postoperative infection rates and lymphatic complications. Functional assays will also be critical to determine whether restoring regulators like RUNX1 can reverse migratory defects without exacerbating inflammation. This complexity may underlie both the heightened inflammatory tone and reduced immune resilience seen in SV/HLHS patients.