Whole blood was collected in PAXgene Blood RNA tubes (PreAnalytiX) and stored at −80 °C prior to extraction. Skin fibroblasts were harvested, and total RNA was isolated using the QIAsymphony system (Qiagen). RNA integrity and concentration were assessed using the Agilent TapeStation, according the manufacturer’s instructions.

Library preparation was performed at the Centre for Genomics Facilities (CFG) under ISO 15189–certified conditions. Polyadenylated RNA was enriched using the KAPA mRNA HyperPrep Kit (Roche) and converted into cDNA libraries following the manufacturer’s instructions. Library quality was evaluated by Qubit fluorometry and fragment-size profiling using an Agilent Bioanalyzer or TapeStation. Paired-end 150 bp sequencing was performed on the Illumina NovaSeq 6,000, targeting approximately 30 million read pairs per sample. Raw reads underwent quality assessment with FastQC and adapter/quality trimming with fastp. Libraries failing QC requirements were excluded from downstream analyses.

Quality-filtered reads were processed with fastp, including adapter removal, low-quality base trimming, and filtering of short fragments. Transcript quantification was performed using Salmon v1.10.1 in quasi-mapping mode against the GRCh38 reference transcriptome and summarized to gene-level counts using tximport.

In the blood whole cohort, sequencing depth averaged ∼104 million 151-bp paired-end reads, with GC content between 49% and 55%. Most libraries exhibited duplication levels consistent with poly(A)-selected RNA. Two samples (Schaefer and Morgan, 2011; Nicholson and Suresh Kumar, 2011) failed RNA-seq quality control and were excluded from downstream analyses based on concordant deviations across multiple technical metrics.

Sample 6 showed markedly elevated duplication rates (86.8%), reduced library complexity, and the lowest number of detected genes in the cohort, while sample 5 exhibited borderline library quality characterized by relatively high duplication rates, and fewer detected genes compared with the remaining samples. Importantly, both samples displayed the HAFOUS episignature, indicating that they do not represent a biologically distinct subgroup.

Principal component analysis (PCA) and sample-sample correlation plot (Supplementary Figure S2) reflected these technical differences and were used as supportive diagnostics rather than as exclusion criteria. Complete RNA and RNA-seq quality control metrics are provided in Supplementary Table 2a.

All four fibroblast case and four control libraries met QC criteria. GC content ranged from 47% to 51%, and sequencing depth ranged between 45 and 90 million reads. High duplication levels were observed, as expected for deeply sequenced fibroblast RNA with limited transcriptomic diversity; however, all samples showed sufficient gene detection and coherent clustering and were therefore retained.

Complete QC metrics are provided in Supplementary Tables 2a,b, with graphical summaries in Supplementary Figures S2, S3.

Gene-level counts were analyzed with DESeq2 (v1.44) (Love et al., 2014). Transcripts reflecting fewer than 10 total counts across all samples were excluded to improve statistical power. Counts were normalized using the median-of-ratios method implemented in DESeq2. For fibroblast analyses, models included age, sex, and age at biopsy as covariates. Differentially expressed genes (DEGs) were defined as those with a Benjamini–Hochberg adjusted p < 0.05. For blood analyses, the same pipeline was applied, with models adjusted for age and sex.

GO Biological Process enrichment was performed separately for blood and fibroblast DEGs using clusterProfiler (Yu, 2024), with gene IDs mapped to Entrez IDs via org. Hs.e.g.,.db. GO with Benjamini–Hochberg adjusted p < 0.05 were considered significant. Gene–concept network plots were generated. Transcription factor (TF) activity was inferred per tissue with decoupleR (Badia et al., 2022) (ULM method) using DoRothEA (Garcia-Alonso et al., 2019) regulons, and visualized via volcano plots. The ULM method estimates TF activity by modeling the relationship between target gene expression and TF–target weights derived from curated regulatory networks. DoRothEA provides confidence-weighted TF–target interactions inferred from literature and large-scale datasets. Resulting TF activity scores and pseudo p-values were visualized as volcano plots.

For the blood interative sub-cohort, four patients and four controls with high-quality matched DNA methylation and RNA-seq data were included in the eQTM analysis. All four cases and four controls from the fibroblast cohort had matched datasets and were similarly analysed. Cis-eQTM analyses were performed separately for each tissue.

Pearson correlation coefficients were calculated between the mean methylation level (β values) of each DMR and the log-transformed expression of its associated gene within a ±250 kb window (DMR-linked genes). This distance was chosen as it reflects the typical range of gene regulatory interactions (e.g., enhancer) established in large-scale resources such as GTEx and Blueprint Epigenome (Ruiz-Arenas et al., 2022; Küpers et al., 2022). Statistical significance was assessed by comparing observed correlations with a null distribution generated from 1,000 permutations of length-matched random DMRs. Bootstrapped 95% confidence intervals were computed, and associations with p < 0.05 were considered significant. The expression of our genes of interest (significant DEGs/eQTMs), across relevant tissues, including brain, skin, skeletal muscle, liver, and whole blood, were evaluated using GTEx v10 median TPM values.

Differentially methylated regions (DMRs) identified from Blood whole cohort and Fibroblast cohorts mapped to hg38 and converted to hg19 using UCSC LiftOver, retaining only uniquely mapped intervals. Public ChIP-Rx datasets for BCOR, H2AK119ub1, and H3K27me3 from the USP7- haploinsufficient neuroblastoma model (GSE263491) were downloaded as bigWig files and imported using rtracklayer. All ranges were harmonized to NCBI chromosome notation. For each replicate, high-signal “peak-like” regions were defined by selecting continuous tiles with ChIP-Rx signal above the 75th percentile; replicates for each mark were merged by concatenation. To test whether DMRs were enriched for Polycomb-associated chromatin features, we generated chromosome- and length-matched background regions using regioneR: randomizeRegions. Background was additionally restricted to CpG islands (UCSC CpGIslandExt) to match CpG density between DMRs and background regions. Average ChIP-Rx signal between DMRs and background were compared using one-sided Wilcoxon rank-sum tests. Overlap was quantified by counting both (i) high-signal tiles overlapping ≥1 DMR and (ii) DMRs overlapping ≥1 tile. Distances between each DMR and the nearest high-signal BCOR, H2AK119ub1, or H3K27me3 region were computed using GenomicRanges: distanceToNearest and log-transformed for visualization. Combinatorial overlap patterns were summarized using UpSet plots. Genomic visualizations were produced using Gviz with EnsDb. Hsapiens. v75 gene annotation alongside DMR and Polycomb tracks.

To assess the clinical and functional relevance of DMR-linked genes identified in blood and fibroblast cohorts, we performed a literature and database search of 37 candidate genes. For each gene, we searched public databases including OMIM (https://omim.org), GeneCards (https://www.genecards.org) and queried PubMed (September 2025) using gene names and functional aliases in combination with neurodevelopmental and brain-related keywords (i.g., “brain”, “neurodevelopmental disorder”, “brain development”, “cerebral cortex”, “hippocampus”, “prefrontal cortex”) as MeSH terms or in free text [].

Our PubMed search strategy was as follows:

(GENE [tiab] OR “GENE full name” [tiab] OR “GENE gene” [tiab]) AND (“Brain” [Mesh] OR brain [tiab] OR “Neurodevelopmental Disorders” [Mesh] OR neurodevelopment*[tiab] OR “brain development” [tiab] OR “Cerebral Cortex” [Mesh] OR “Hippocampus” [Mesh] OR “Prefrontal Cortex” [Mesh]).

Where appropriate, gene-specific functional keywords or disease associations were added, based on information from OMIM or GeneCards. For each gene, clinical phenotypes and functional roles reported in the literature were documented, and the presence or absence of a known Mendelian phenotype was annotated. The resulting data were summarized in a table for further interpretation.

The blood whole cohort comprised nine patients carrying pathogenic or likely pathogenic USP7 variants and four healthy controls (Table 1). All nine patients exhibited the characteristic HAFOUS methylation episignature in peripheral blood, including six patients previously reported in our original HAFOUS episignature study. High-quality DNA methylation data were obtained for all participants, forming the blood whole cohort used for differential methylation and DMR analyses.

RNA sequencing was performed for six out of the nine USP7 cases together with the same four healthy controls as used for the DNA methylation. Two case samples showed insufficient library complexity and elevated duplication rates and were excluded following RNA-seq quality control. The remaining four USP7 cases and four healthy controls constituted the blood integrative sub-cohort, which was used for differential gene expression and eQTM analyses.

Comprehensive demographic, molecular, and analytic inclusion details for all samples are provided in Table 1.

Genome-wide analysis identified 209 significant differentially methylated positions at FDR <0.05, of which 202 (96.7%) were hypermethylated and 7 (3.3%) were hypomethylated. When examining the top 5,000 DMPs (ranked by adjusted P-value), 4,379 (87.5%) were hypermethylated and 621 (12.5%) were hypomethylated, confirming a pronounced global hypermethylation pattern in the cohort. The top 5,000 DMPs ranked by effect size are provided in Supplementary Table 3a.

DMRcate identified 17 significant DMRs (≥5 CpGs, 1 kb window, FDR <0.05), all showing consistent hypermethylation (Supplementary Table 4a). To ensure that methylation patterns were comparable between the full cohort and the integrative sub-cohort, we compared mean β-value distributions for each DMR across both groups. Wilcoxon rank-sum tests showed no meaningful differences between whole-cohort and sub-cohort methylation distributions (all p > 0.05 except one control comparison), with uniformly small effect sizes (Cliff’s delta), confirming that the sub-cohort is representative of the whole cohort (Supplementary Table 4d).

After filtering out lowly expressed genes (fewer than 10 total counts across all samples) and retaining only genes expressed in at least 75% of both case and control samples, 15,295 genes remained for differential expression analysis. RNA-seq analysis of whole blood (n = 4 cases; n = 4 controls) identified 33 differentially expressed genes (DEGs) meeting the Benjamini–Hochberg adjusted p < 0.05 significance threshold (Supplementary Table 5a). Only 33 genes were significant and Log₂ fold changes ranged from approximately −4 to −30, primarily reflecting genes with low or absent expression in one group, as expected given the small cohort size.

Gene set enrichment analysis (GSEA) performed on the ranked blood transcriptome did not identify any gene ontology categories or pathways reaching significance after multiple testing correction. The lack of enrichment is consistent with the limited number of differentially expressed genes and the overall modest transcriptional alterations observed in blood.

Using the 17 blood-derived DMRs, cis-eQTM analysis within a ±250 kb window identified nine significant DMR–gene associations (p < 0.05) (Supplementary Table 6a). These pairs showed consistent methylation–expression relationships indicative of local regulatory effects. Four associations demonstrated strong correlations (r ≥ 0.7) and are presented in Figure 1, whereas the remaining significant associations with moderate correlations (r < 0.7) are shown in Supplementary Figure S3.

Primary skin fibroblasts were obtained from four patients with (likely) pathogenic USP7 variants or deletions and four unaffected controls (Table 2). All samples generated high-quality DNA methylation and RNA-seq data and were included in downstream analyses. The same eight samples therefore formed both the fibroblast cohort for methylation analysis and the fibroblast integrative cohort for transcriptome and eQTM analyses.

Differential methylation analysis of fibroblast samples identified no DMPs that reached FDR significance after multiple-testing correction, consistent with the small cohort size. Therefore, we proceeded with exploratory interpretation based on nominal p-values. The top 5,000 CpG sites ranked by nominal significance are summarized in Supplementary Table 3b. Of these, 2,507 (50.1%) were hypermethylated and 2,493 (49.9%) hypomethylated in USP7-deficient fibroblasts compared with controls, indicating an approximately balanced global distribution of methylation changes. These DMPs represent candidate loci for downstream region-based and integrative analyses, including DMR identification and eQTM correlation.

Using bumphunter (cutoff = 0.15; 100 bootstraps), we identified 2,143 nominally significant DMRs in fibroblasts, comprising 1,314 (61.3%) hypermethylated and 829 (38.7%) hypomethylated regions. No regions passed FWER correction, which is expected given the cohort size and the bumphunter framework, but the set of nominal DMRs was stable across bootstrap iterations. These DMRs formed the basis for downstream integrative analyses, including eQTMs and chromatin contextualization.

To investigate transcriptional alterations associated with HAFOUS, we performed differential gene expression analysis comparing fibroblasts derived from affected patients to healthy controls. We observed a total of 310 genes that were significantly differentially expressed (BH adjusted p < 0.05, |log₂FC| >0), comprising 165 upregulated and 145 downregulated genes in patient samples. The expression profiles of these genes with full results are provided in Supplementary Table 5b.

Unbiased gene ontology enrichment analysis of differentially expressed genes in USP7-patient fibroblasts revealed significant enrichment in biological processes related to developmental and cellular pathways, notably those associated with nervous system development, including positive regulation of cell development, axon guidance, regionalization, and peptide transport. The gene–concept network visualization (Figure 2A) highlighted key genes contributing to multiple of these processes, such as FGF2, RELN and BMP2.

To investigate upstream regulatory influences, transcription factor activity was inferred using the decoupleR ULM method with DoRothEA regulons. Several TFs showed significantly altered activity scores (Figure 2B), including increased activity of CREB1 and RELA, and decreased activity of ASCL1, E2F1, and KLF5. These activity changes reflect coordinated shifts in target gene expression rather than changes in TF transcript levels.

As stated earlier, we identified 2,143 differentially methylated regions (DMRs) between patient fibroblasts and controls using the bumphunter algorithm (Δβ cutoff = 0.15; 100 bootstrap iterations). To investigate whether USP7 haploinsufficiency disrupts transcription through epigenetic mechanisms, we performed an integrative expression quantitative trait methylation (eQTM) analysis, linking DMRs to gene expression within a ±250 kb cis-window. This analysis uncovered 559 significant eQTMs (Supplementary Table 6b). Gene set enrichment of associated loci revealed Biological processes related to developmental and morphogenetic processes, particularly anterior–posterior patterning and skeletal system development, driven largely by HOX gene clusters and embryonic developmental regulators (Supplementary Table 7a). Among these, 16 genes were also differentially expressed, pinpointing candidates under transcriptional and epigenetic regulation (Supplementary Table 7b). Expression analysis confirmed marked downregulation of HOXB3, HOXB5, and HOXB6 in patient fibroblasts compared with controls (Figure 3A). Quantification of DMRs per gene revealed that the HOX cluster, together with CHRD, carried the highest DMR burden, underscoring their potential role as regulatory hotspots (Figure 3B).

A global integration plot demonstrated two distinct regulatory patterns among the 16 overlapping eQTMs (Figure 3C):

HOX cluster genes (HOXB3, HOXB5, HOXB6) showed strong negative correlations (r < −0.7, log₂FC < −2), consistent with hypermethylation-associated repression.

Within the HOXB cluster, HOXB3 is shown as an example. Six DMRs within ±250 kb of HOXB3 were significantly correlated with expression, including five strong negative correlations (r = −0.76 to −0.87, p < 0.05) consistent with hypermethylation-mediated silencing, and one positive correlation (r = 0.83, p = 0.015) (Figure 3D). All significant DMR–gene associations and their corresponding correlation statistics are provided in Supplementary Table 6b, and patient DMR–expression correlation plots are presented in Supplementary Figures S5, S6 for visual assessment of association directionality and strength. Baseline tissue-expression profiles from the GTEx database showed that several genes dysregulated in USP7-haploinsufficient fibroblasts including CHRD, SLC7A5, and PSMD2 display moderate to high expression levels in skin and multiple brain tissues (Supplementary Figure S7), supporting their potential relevance to both peripheral and neurodevelopmental features of USP7-associated syndromes.

To assess whether USP7-associated methylation changes are conserved across tissues, we compared fibroblast DMRs (bumphunter; n = 2,143) with blood-derived DMRs identified by DMRcate (n = 17). This analysis identified eight exact overlapping DMRs, located on chromosomes 2, 3, 6, 7, 10, 11, and 16 (Supplementary Table 7c). No additional matches were detected when expanding the window by ±250 bp (Figure 4).

At the gene level, annotation of DMRs within a ±250 kb window revealed five shared genes: NKAPL, GPR75/ASB3, LRRC15, ZNF423, and LRMDA. These reflect cases where distinct DMRs in each tissue map near the same gene, rather than representing shared genomic regions.

Integration with expression demonstrated limited cross-tissue concordance. Fibroblasts yielded 559 significant eQTMs, whereas blood yielded nine, and only two genes (HOXA9 and CNTD1) were common to both datasets. No differentially expressed genes were shared between tissues.

Together, these findings indicate that only a small fraction of USP7-associated methylation changes are conserved across blood and fibroblasts, and that methylation–expression coupling appears largely tissue-specific.

To assess whether differential methylation in USP7-haploinsufficient fibroblasts preferentially arises within Polycomb-enriched chromatin, we intersected all 2,143 DMRs with BCOR, H2AK119ub1, and H3K27me3 ChIP-Rx profiles generated in a USP7- haploinsufficient neuroblastoma model (Wolf van der Meer et al., 2025). For both hyper- and hypomethylated DMRs, ChIP-Rx signal was significantly higher than matched, chromosome, length-controlled and cpg background regions, indicating robust enrichment of Polycomb-associated features (Wilcoxon rank-sum tests; Supplementary Table 8a). Among hypermethylated DMRs, enrichment was strongest for H2AK119ub1 (W = 1,275,865.5; p = 3.6 × 10⁻¹⁰⁰) and BCOR (W = 1,046,247; p = 2.6 × 10⁻²¹), with H3K27me3 showing a weaker but still significant increase (W = 960,827; p = 2.7 × 10⁻⁷). Hypomethylated DMRs exhibited similar patterns (H2AK119ub1: W = 521,745; BCOR: W = 432,786; H3K27me3: W = 378,845; all p < 10⁻³).

Combinatorial overlap analysis (Figure 5) further demonstrated that many DMRs co-localized with multiple Polycomb features. Among hypomethylated DMRs (n = 829; Figure 5A), 165 overlapped all three marks, 127 overlapped H3K27me3 + H2AK119ub1, and 139 uniquely overlapped H2AK119ub1. Hypermethylated DMRs (n = 1,314; Figure 5B) showed comparable patterns, including 251 overlapping H2AK119ub1 + BCOR, 246 overlapping H2AK119ub1 + H3K27me3, and 105 intersecting all three marks (Supplementary Tables 8c–e).

Genome-wide quantification supported these observations: 72% of all DMRs overlapped high-signal BCOR, 77% overlapped H2AK119ub1, and 26% overlapped H3K27me3 (Supplementary Table 8b). In contrast, only 0.07%–0.14% of all genome-wide high-signal Polycomb tiles overlapped a DMR, indicating that DMRs occupy a highly restricted subset of Polycomb-rich domains.

To assess spatial proximity independent of direct overlap, we measured the distance from each DMR to the nearest high-signal Polycomb region. Both hyper- and hypomethylated DMRs were located within tens to hundreds of base pairs from PRC1-associated BCOR and H2AK119ub1 peaks, whereas distances to H3K27me3 were broader (Figure 6), consistent with the diffuse nature of PRC2 occupancy.

Finally, applying the same analysis to the 17 hypermethylated DMRs identified in blood showed a qualitatively similar but more modest pattern. Blood DMRs were also enriched for BCOR and H2AK119ub1 signal relative to background, whereas H3K27me3 showed no significant increase (Supplementary Table 9a). Consistent with fibroblasts, most blood DMRs overlapped PRC1-associated regions (BCOR: 82%; H2AK119ub1: 76%), while a minority overlapped H3K27me3 (18%) (Supplementary Tables 9b–e).

To contextualize the regulatory loci identified in our methylation and eQTM analyses, we reviewed clinical and functional annotations for 37 genes associated with HAFOUS-linked DMRs using OMIM and the published literature (Table 3). A minority of these genes have established Mendelian phenotypes; examples include UNG (immunodeficiency with hyper-IgM type 5), HAGH (glyoxalase II deficiency), and GPR68 (amelogenesis imperfecta). Additional entries such as MARCHF6, reported in familial adult myoclonic epilepsy, and DAPK1, a kinase involved in neuronal apoptosis and plasticity, further highlight links between the identified loci and known disease or neuronal pathways.

For most genes, no Mendelian disorder has been described, yet their biological functions point to plausible regulatory relevance. These include early developmental transcription factors (HOXA9, HOXB3, HOXB6, ZKSCAN4, ZNF165), components of ribosome biogenesis (BRIX1, LTO1), DNA repair (NEIL2), and signaling pathways such as WNT and Hippo (CCN4, WWC2). Additional genes with roles in neurodevelopmental processes such as SLC7A5, an amino acid transporter linked to autism spectrum traits, and CHRD, a BMP antagonist essential for early forebrain induction provide further functional context. Among these, MARCHF6 stands out given its involvement in ubiquitination pathways that intersect conceptually with USP7 function, while the presence of multiple HOX genes (HOXB3, HOXB5, HOXB6) aligns with the broader domain-level regulatory disruptions identified in fibroblasts.