2. Methods

A brief overview of methods is provided below; full details for all analyses are available in File S1.

2.1. Patient Recruitment and Ethics

Study participants were recruited at Cincinnati Children′s Hospital Medical Center, Ohio, United States. A total of 14 trios and one duo were enrolled, comprising probands with a confirmed diagnosis of NDD and/or CHD and their biological parents. Participants were consecutively enrolled without restriction by phenotype, severity or genetic aetiology, forming a heterogeneous, clinically referred pilot cohort. Written informed consent was obtained from all participants, including parental consent for minors. The study was approved by the Cincinnati Children′s Hospital Institutional Review Board (IRB ID 2021‐0493).

2.2. Sample Collection and Sequencing

Peripheral blood samples were collected from all participants. Genomic DNA and total RNA were extracted using standardised protocols. WGS was performed on 13 trios and one duo at ≥ 30x coverage using Illumina sequencing at the Broad Institute of MIT and Harvard. RNA sequencing was performed on a subset of nine trios and one duo (≥ 25 million paired‐end reads), and DNA methylation profiling was performed on 14 trios and one duo using the Illumina Infinium MethylationEPIC BeadChip array. Standard WGS quality control metrics (relatedness, sex confirmation, coverage, contamination and insert size) were assessed; all samples passed. The availability of whole‐genome sequencing, RNA sequencing and DNA methylation data for each participant is summarised in Table S1.

2.3. Annotation of Sequencing Data

WGS data were preprocessed and analysed according to the GATK Best Practices workflow [10]. Reads were aligned to GRCh38 using BWA‐MEM (v0.7.17) [11], variants called using GATK HaplotypeCaller with joint genotyping and quality filtered using variant quality score recalibration. Variants were annotated using ANNOVAR [12] and the Ensembl Variant Effect Predictor [13].

2.4. Variant Prioritisation

Variants with MAF < 1_%_ in gnomAD (v2 and v3) [14] were assessed. Variants were considered potentially damaging if they altered protein length (nonsense, frameshift, in‐frame indel or splice site), were missense variants with REVEL ≥ 0.6 [15] or BayesDel ≥ 0.1 [16] or were exonic structural variants. Predicted‐damaging variants in NDD/CHD–associated genes (Table S2) were evaluated and verified by visual inspection using IGV [17]. For CHD trios, coding variants across all genes were additionally assessed under autosomal dominant (MAF < 0.001), autosomal recessive (MAF < 0.01), de novo and compound heterozygous inheritance models, excluding highly polymorphic genes. Variants of interest were classified according to ACMG guidelines. Genome‐wide high‐confidence LoF variants were identified using LOFTEE [14], filtered to no LoF flags, allele count < 20 in gnomAD 4.1, present in < 5 cases, read depth ≥ 8, allelic balance ≥ 0.15 and annotated to LoF‐intolerant genes (LOEUF ≤ 0.35).

2.5. Polygenic Risk Score Analysis

PRS were calculated for the nine probands of European ancestry using the pgsc_calc reproducible workflow [18], relative to a healthy elderly European control cohort (n = 3823) [19]. Scores were computed for ADHD (PGS003753) [20], anxiety (PGS001021) [21], autism (PGS000327) [22] and depression (PGS001829) [23]. PRS were transformed to control distribution percentiles; group differences were assessed using a Mann–Whitney U test with Benjamini–Hochberg correction.

2.6. RNA Sequencing Analysis

Reads were aligned to GRCh38 using STAR (v2.7.10) with two‐pass alignment. Aberrant expression was identified using OUTRIDER (v1.8.2) [24], supplemented with 97 RNA sequencing samples from children with complex medical disorders to improve outlier detection [25]. Six samples with size factors outside the 0.7–1.2 required range were excluded; the remaining 24 samples were analysed. Genes with adjusted p value < 0.05 were flagged as aberrantly expressed. Splicing variant validation and mono‐allelic expression (MAE) analysis were performed on the 10 trios with available transcriptomic data.

2.7. DNA Methylation Analysis

Methylation data were quality‐controlled using standard criteria; probes failing detection p value, bead count, cross‐reactivity, CpG‐SNP or sex chromosome filters were removed. One proband failed QC and was excluded. Beta values were normalised using functional normalisation. Epigenetic age was estimated using the Hannum blood‐specific and Horvath pan‐tissue clocks, which are optimised for blood samples and paediatric/adult cohorts, respectively; age acceleration was defined as the difference between predicted and chronological age. To account for family structure, linear mixed‐effects models were fitted with family as a random effect and age and sex as covariates; group differences were assessed using ANOVA with pairwise t‐tests.

2.8. Cross‐Cohort Differential DNA Methylation Analysis

To enable individual patient‐level methylation comparisons, the study cohort was integrated with two GEO reference control datasets: GSE74432 (n = 53, Illumina 450K [26] and GSE97362 (n = 125, Illumina 450K) [27]. Raw IDAT files were processed using minfi (v1.46.0) with preprocessNoob normalisation and dye bias correction. Batch effects were corrected using ComBat (sva v3.48.0) with negative control probe intensities as covariates; effectiveness was confirmed by PCA (Figure S1). Cell type composition was estimated using EpiDISH (v2.16.0) with the IDOL whole‐blood reference panel. A secondary epigenetic age analysis was performed using the EN and Levine clocks (selected for CpG coverage), with group differences assessed using ANOVA and Tukey′s HSD. For each proband, differentially methylated probes (DMPs) were identified using limma (v3.56.2) with design matrices incorporating group status, batch and cell type proportions, applying empirical Bayes moderation (FDR < 0.05). Probands with > 10,000 significant DMPs were classified as having methylation dysregulation.

2.9. Methylation Pathway Analysis

Within‐cohort differentially methylated regions (DMRs) were identified using the MethylMiner pipeline [28], (≥ 5 probes within a 500‐bp window). Pathway enrichment for genes annotated to proband‐specific DMRs was assessed using clusterProfiler [29] with Gene Ontology Biological Process terms [30, 31].