All eight centers for human genetics in Belgium engaged in this prospective RCT – the BeSolveRD trial – with the primary endpoint of comparing the use of SoC and GS for the molecular diagnosis of DD including ID and CA. The study was registered (ClinicalTrials.gov NCT07051213), approved by the Ethical Commission of UZ Leuven (S64603) and was performed in accordance with the Declaration of Helsinki. Informed consent was obtained from all participants or their legal representatives. Inclusion criteria were (1) moderate to profound ID, (2) mild to moderate ID with recurrence among siblings and parents with a normal intellect, (3) mild to moderate ID and dysmorphism, (4) mild to profound ID and one or more CA, (5) one major CA and dysmorphism and, (6) multiple major CA in two or more different organ systems. CA are considered major when they have a significant impact on the health of the child, whereas rare morphological features with a limited functional or esthetic impact are considered as minor anomalies [26]. Dysmorphism was defined as the presence of at least 3 minor anomalies. Participants were included without age-related restrictions. Probands with a suspicion of an acquired cause, e.g. congenital infection and toxic exposure, were excluded. Probands with prior testing of a gene panel targeting multiple conditions or exome analyses were also excluded. Prior CMA or sGS was not an exclusion criterion. Participants were recruited by clinical geneticists to ensure detailed phenotypic assessments were performed and thorough family histories were collected, both being essential for variant interpretation.

Blood was collected from probands and their parents to enable trio analyses with inheritance based-variant filtering. For 16 probands only one parent was available, no singletons were included. Nineteen pairs of siblings with a similar DD and their parents were also included. The initial aim was to recruit 800 samples. The recruitment process was slower than anticipated due to an overestimation of both the proportion of patients meeting the inclusion criteria and the study’s participation rate. The recruitment was ended after 3 years (02/06/2021–05/07/2024) with a total of 567 participants included. The lower sample size impacts the power of the study but does not introduce a bias between both study arms given the prospective and randomized repartition of included samples, as well as the recruitment stop determined before completion of the data analysis. Participants were randomized using a blocked randomization list [27] into the SoC arm of the study (n = 283), or the GS arm (n = 284), with center-specific randomization lists used to ensure balanced allocation across study arms within each laboratory (Additional File 1: Table S1). The randomization resulted in a similar distribution of demographic characteristics and indications categories between both study arms (Fig. 1). Data collection and study management were conducted using REDCap [28].

In the SoC arm, all probands were investigated using the standard methodology employed at each participating center: CMA- (5/8 centers) or sGS- (3/8 centers) based CNV analysis and subsequent ES (Additional File 1: Table S2). In the different centers, the CMA resolution varied between 50 kb and 200 kb with higher resolution in genes of interest and for sGS the coverage varied between 1x and 1.5x. ES was performed following the wet lab and bioinformatic protocols that had been previously validated for clinical use by each of the 8 centers [29]. The wet lab protocols mainly differed in the shearing methodologies (enzymatic (n = 5), sonication (n = 3)), and the exome capture kit used (KAPA HyperCap (n = 2), Twist exome (n = 5), xGen exome (n = 1)). Differences in bioinformatic processing included the version of the reference genome used (hg19 vs. hg38 in three and five centers, respectively), pre-processing of the reads prior to mapping (performed in two centers), and base quality score recalibration (conducted in five centers). The tools used for mapping and duplicate marking also differed.

For the GS arm, library preparation and sequencing were performed in five centers for the eight recruiting centers [29]. DNA was sheared using different methods, depending on the site, including sonication (n = 2), tagmentation (n = 2) and enzymatic fragmentation (n = 1). Libraries were prepared without PCR amplification using three distinct kits: Roche KAPA HyperPrep PCR-free (n = 2), Illumina DNA PCR-free (n = 2), or NEBNext Ultra II FS DNA PCR-free (n = 1). Paired-end sequencing was performed on Illumina NovaSeq6000 platforms. Bioinformatic processing was done by each of the eight recruiting centers using a custom pipeline consisting of mapping to the reference genome, duplicate marking and variant calling. Major differences in the analysis pipelines included the version of the reference genome (hg19 vs. hg38 in two and six centers, respectively), the pre-processing of the reads prior to mapping (done in one center) and the base quality score recalibration (done in five centers). The tools used for duplicate marking and base quality score recalibration also differed. Six of eight centers extracted an exome VCF file prior to interpretation. All centers but two called CNVs from GS data. A mean informative coverage of 42X (+/- 10%) computed by either Picard [30] or mosdepth [31] was targeted. Uniformized re-analysis of the data of GS negative samples, including non-coding regions, was performed across the eight centers with the Emedgene software v.32.0.29 (Illumina Inc., San Diego, CA, USA) [32]. Beyond SNVs/indels, the software allows the detection of short tandem repeats (STRs) and SVs.

SoC methods were validated prior to the study. To ensure a swift implementation and validation of GS in all participating centers, all centers engaged in an inter-laboratory ring trial before patient recruitment. Two Genome-In-A-Bottle (GIAB) samples [33], HG001 and HG002, were sequenced in the five sequencing centers. Data was shared with the eight analyzing centers so that each center could analyze the samples with its own pipeline and compute coverage statistics, precision, and recall with a standardized protocol. Comparison of the results helped each center pinpoint problematic steps in GS library preparation and sequencing, highlighting the need for large insert sizes and optimal pooling of the libraries to avoid optical duplicates and obtain the expected mean sequencing depth of all samples. The use of different bioinformatic analysis pipelines resulted in only minor differences in the performance of SNV and indel calling as similar precision (HG001: 0.9669 (+/- 0.0229); HG002 0.9873 (+/-0.0076)) and recall (HG001: 0.9841 (+/- 0.0067); HG002 0.9898 (+/-0.0074)) were obtained in all centers (Additional File 1: Fig S1).

In the SoC arm, variant interpretation was performed within the accredited diagnostic laboratories of each participating center using their established software platforms and filtering pipelines (Additional File 1: Table S2). Identified variants were prioritized based on a gene panel; in three of the eight centers additional exome-wide filtering based on the phenotype of the patient using Human Phenotype Ontology (HPO) terms was applied. In all centers, the variants were further assessed according to Mendelian inheritance modes, starting with de novo variants and further investigating X-linked (XL) variants as well as variants that might lead to autosomal recessive (AR) conditions (homozygous and compound heterozygous variants). In addition to these standard analyses, potential disease-causing inherited variants in genes associated with autosomal dominant (AD) conditions were also investigated. Variants were classified according to the American College of Medical Genetics and Genomics/Association for Molecular Pathology and/or Association for Clinical Genomic Science classification guidelines [34, 35]. (Likely) pathogenic variants (classes 4 and 5) that could at least partially explain the phenotype of the patient were reported and considered as certain diagnoses. VUS were only reported upon high suspicion of association to the patient’s phenotype. Such VUS, as well as (likely) pathogenic variants with an uncertain contribution to the phenotype, mostly susceptibility factors, were considered as uncertain diagnoses.

In the GS arm of the study, variant filtering and interpretation were performed by the investigators of the same diagnostic laboratories to avoid interpretation biases. In the first phase, GS SNVs/indels as well as CNVs (for six centers), were interpreted using the same software, filters and interpretation criteria as in the SoC arm. Additionally, re-analysis of the data of negative GS cases was performed using the Emedgene software [32] to tackle potential causal repeat expansions, structural and/or intronic variants. The AI-based filtering and interpretation of variants in this software integrates HPO-based phenotypic information to improve variant prioritization [36].

Comparative analyses of diagnostic yields were conducted upon completion of recruitment and finalization of variant interpretation in the different diagnostic laboratories. The observed difference in diagnostic yield between the study arms was adjusted to account for potential biases, including the inclusion of siblings and differences in variant filtering strategies and analysis timing resulting from the Emedgene reanalysis. The impact of additional confounding factors, such as cohort demographic characteristics, was also assessed. Adjustment for sex distribution was performed using logistic regression. Statistical significance of the differences in diagnostic yield was measured using the score method with skewness correction, based on the approaches of Laud, Gart and Nam [37, 38].

A total of 567 probands were recruited into the study, including 353 males and 214 females. The randomization process resulted in a similar distribution of demographic characteristics and clinical indications between the study arms (Fig. 1). Males accounted for 65.7% (186/283) of the SoC arm and 59.2% (168/284) of the GS arm. Most probands were either toddlers (1–3 years old; 156/567, 27.5%) or children (4–11 years old; 207/567, 36.5%), while smaller proportions were infants between one month and one year of age (28/567, 4.9%), adolescents (12–17 years old; 75/567, 13.2%), or adults (≥ 18 years; 93/567, 16.4%). Given the low uptake of the study in neonates (< 1 months old; 3/567, 0.5%) and in the prenatal setting (products of conception; 5/567, 0.9%), the study findings primarily reflect diagnostic yields in probands older than one month. Moderate to profound ID (268/567, 47.3%) and mild to moderate ID in association with dysmorphism (177/567, 31.2%) constituted the two most frequent indications for genetic testing. The four other indications each contributed for 4.4% to 6.7% (25–38 probands) of the overall cohort. For 377/567 probands a CNV analysis had been performed in the past without identifying an underlying causal variant (204/283, 72.1% in the SoC arm and 173/284, 60.9% in the GS arm). Fragile-X syndrome was tested before inclusion for 159/353 males and 65/214 females and a conventional karyotype had been performed for 90 probands. Sanger sequencing had previously been performed to assess a single gene or a disease-targeted panel of genes for respectively 39 and eight individuals, and NGS had been performed to evaluate a single gene or a disease-targeted panel of genes in 19 and 61 individuals, respectively (Fig. 1).

A certain molecular diagnosis was established for 76/283 (26.9%) probands analyzed with SoC methodologies (evaluating de novo, XL and AR variants). Most of these variants (70/76, 92.1%) were SNVs or indels identified by ES; the remainder were identified using CMA/sGS and included one SV, four SNVs and one aneuploidy (45,X). Uncertain diagnoses were identified in 26/283 (9.1%) of the probands. In the GS arm (evaluating de novo, XL, and AR variants), a certain molecular diagnosis was made for 100/284 probands (35.2%). Similarly as in the SoC arm, the majority (92/100, 92.0%) were SNVs or indels. Seven variants were CNVs and one STR expansion was identified. In addition, uncertain diagnoses (variants of unknown significance and/or unknown contribution to the phenotype) were identified for 27/284 (9.5%) of probands (Fig. 2, Additional File 1: Fig S2, Additional File 2: Table S3).

Eight of the certain diagnoses in the GS arm were detected during the systematic reanalysis using the Emedgene platform. One of these was a de novo 5’UTR SNV in MEF2C (MIM: 600662) NM_002397.5:c.-26 C > T. This non-coding variant has been shown to create an upstream new start codon which leads to the translation of an elongated protein [8]. In another proband, a STR expansion in the promotor region of FMR1 (MIM: 309550), responsible for Fragile-X syndrome, was identified. In two other, unrelated probands, a pathogenic NM_033517.1:c.3679dupG p.(Ala1227Glyfs*69) in SHANK3 (MIM: 606230) had been filtered out in the routine analysis because of its frequency in the GnomAD and Exome Sequencing Project population databases, which might be inflated due to sequencing errors as the insertion is affecting a stretch of eight guanines [39]. A fifth variant was NM_145686.4:c.572C > A p.(Thr191Asn) in MAP4K4 (MIM: 604666). This variant was initially classified as a VUS but reclassified as pathogenic during Emedgene reanalysis due to novel evidence [40]. Three additional diagnoses were uncovered because the applied gene panel did not include the involved genes (INTS11 (MIM: 611354) NM_0178716:c.1003C > T (p.Gln335Ter), c.1313C > T (p.Pro438Leu); NKX6-2 (MIM: 605955) NM_177400:c.134delTinsGG (p.Leu45ArgfsTer415); H3F3A (MIM: 601128) NM_002107:c.365C > T (p.Pro122Leu)).

Beyond the standard analysis of de novo, XL, and AR variants, an additional assessment was performed to identify inherited variants in genes associated with AD conditions. This analysis identified nine inherited (likely) pathogenic variants in the SoC and 13 in the GS arm. These 22 additional diagnoses represent 11% of all certain diagnoses and an increased diagnostic yield of 3.9% (Figs. 2B and 3). For most of these families (18/22) there was no prior suspicion of a dominantly inherited disorder. Seventeen of the 22 variants (77.3%) were interpreted as fully explaining the patient’s phenotype while five were recurrent CNVs considered risk factors for ID and were thought to partially explain the clinical presentation of the patient. As a comparison, the de novo analyses enabled identification of a (likely) pathogenic variant for 134 patients (23.6%) while the XL analysis was contributive for only 11 patients (1.9%) and the recessive analysis for 23 patients (4.1%). AD inherited variants of unknown significance and/or unknown contribution to the phenotype were retained in another 13 patients (2.3%) (Additional File 1: Fig S3). One of these variants was a ~ 30 kb deletion that includes the 5’UTR first exon of ANKRD11 (MIM: 611192). Deletions of this region have been associated with downregulation of ANKRD11 transcription, resulting in KBG syndrome [41].

The variant allele frequency (VAF) of the inherited SNVs (0.30 to 0.57) was consistent with heterozygous carriership for most (20/22) carrier parents (Additional File 1: Table S4). Two causal variants were likely inherited from a mosaic parent (VAF < 0.25), which is a plausible explanation for the absence of a phenotype in the parent. Another cause of differences in parental and patient phenotypes is imprinting, which are involved in the variant disrupting the GNAS (MIM: 139320) gene [42]. Overall, inheritance was maternal for 15/22 cases (68.2%) and paternal for 7/22 (31.8%). Half of the carrier parents (11/22) were reported to be asymptomatic. In eight, mild symptoms of the disease could be identified, and three presented more severe symptoms. For all genes in which (likely) pathogenic inherited variants were identified, variable expression and reduced penetrance has been described in the literature.

Comparing the diagnostic yield in both arms of the study (evaluating de novo, XL and AR variants), we observed a significant increased yield of 8.4% (95% CI: 0.7% − 15.9%, p = 0.032) in the GS arm (SoC 26.9% (76/283), GS 35.2% (100/284)). After analysis of inherited AD variants, the added value of GS increased to 9.8% (95% CI: 1.9% − 17.5%, p = 0.015) (Table 1), which corresponds to a relative added diagnostic yield of 32.7%.

The impact of SVs/aneuploidies on the yield difference was minor (+ 1.0%), with the number of (likely) pathogenic SVs being low in both arms (seven CNVs and one SV in the SoC arm and 11 CNVs in the GS arm). The added diagnostic yield in the GS arm primarily reflects an increase in the detection of SNVs/indels, from 26.9% in the SoC arm to 35.6% in the GS arm (+ 8.7%, 95% CI: 0.011–0.163, p = 0.025). Of the SNVs/indels detected in the GS arm, only one variant was not disrupting the coding genome or splice sites, namely the 5’UTR NM_002397.5:c.-26C > T variant in MEF2C (MIM: 600662). We investigated whether some of the SNVs/indels identified with GS were located in exonic regions poorly covered in the ES data, but the mean informative coverage of the genomic loci of all GS-detected SNV/indels in the exome data was > 20x (Additional File 3: Table S5). When comparing the yield for each variant subtype, we observed a significant difference in the number of out-of-frame indels, with 26 frameshift variants identified in the GS arm in contrast to 10 in the SoC arm (Fig. 4). This difference corresponds to an additional 5.6% (95% CI: 1.7% − 9,8%, p = 0.006) of patients for whom an out-of-frame indel was identified in the GS arm. GS also enabled diagnosis of one proband with Fragile-X syndrome, one of the most prevalent causes of ID in males [43], which is caused by a tandem repeat expansion in the FMR1 (MIM: 309550) promotor region. Among the CNVs, GS identified a small 5’UTR deletion potentially downregulating ANKRD11 (MIM: 611192), a VUS which would have been missed using SoC.

We further assessed whether certain analytical strategies or characteristics of the cohorts in each of the arms could have impacted diagnostic yields. We reviewed the distribution of the 19 sibling pairs across both study arms and evaluated its potential impact on the results. Excluding one sibling from each pair when both had the same diagnosis and were randomized to the same study arm had only a minor effect on overall outcomes, with diagnostic yields of 30.0% (82/273) in the SoC arm and 40.3% (112/278) in the GS arm (diagnostic yield difference of 10.3%, p = 0.012). We also investigated the potential influence of expanded gene panels, less stringent population-based filtering, and recent literature evidence in the GS Emedgene reanalysis. Correcting for the six molecular diagnoses for which detection was explained by these parameters, we still measured an increase in diagnostic yield of 8.1% (95% CI: 0.2% − 15.9%, p = 0.045) in the GS arm. Given the almost equal distribution of individuals across study arms in each center, the small differences in overall yield between centers are not expected to significantly impact the results of our comparison between both study arms (Additional File 1: Table S1).

We then explored whether differences in yield could be influenced by indication categories or age groups distribution. We mainly observed lower yields for the patients for whom CA are reported without associated ID (25% vs. an average of 36% for the categories including ID). Comparing the yield across the age groups, the yield was lower in children (4–11 years old) than in other age groups. The only difference reaching significance among indication categories or age groups, was the difference in yield between toddlers (65/156, 41.7%) and children (56/206, 27.2%) (difference of 14.5%, 95% CI: 4.8% − 24;4%, p = 0.004). However, distribution of all subcategories was relatively homogeneous amongst study arms (Fig. 1). Comparing the diagnostic yield between males and females, 45.5% (97/213) of females obtained a diagnosis vs. 28.5% (101/354) of males (Figs. 4, Additional File 1: Fig S4, Tables S6&S7), corresponding to a significant difference of 17% (95% CI: 8.7% − 31.7%, p < 0.001). We therefore evaluated to what extent the slight difference in male/female ratio in both study arms might affect our results. After adjusting for this difference, the diagnostic yield gap was estimated at 9.0% (95% CI: 1.3% − 16.6%, p = 0.022) (Additional File 1: Table S6). Combining the potential influence of the difference in yield between both sexes with the corrections for the Emedgene reanalysis reduces the diagnostic yield gap between GS and SoC to 7.3% (95% CI: -0.6% to 15.1%, p = 0.069).

Finally, we evaluated the potential impact of the distribution of autism, a phenotypic trait associated with lower diagnostic yields [44, 45]. In our cohort, the diagnostic yield in probands with reported autism or autistic behavior was 24.1%. Autism or autistic behavior was more frequent among the unsolved (131/306, 42.8%) than the solved cases (45/198, 22.7%) and in males (137/353, 38.8%) than females (54/214, 25.2%), which might explain part of the observed difference in yield between genetic sexes. However, the proportion of probands with autism or autistic behavior was higher in the GS arm (105/284, 37.0%) than the SoC arm (86/283, 30.4%) (Additional File 1: Fig S5). Hence, correcting for this feature would increase the diagnostic yield of GS. Among the other frequent phenotypic traits, higher diagnostic yields were observed in probands with hypotonia (23/39, 59.0%), hypertelorism (16/36, 44.4%), microcephaly (24/55, 43.6%), or short stature (19/44, 43.2%) (Additional File 1: Fig S6). However, no significant differences were observed in the distribution of the frequent phenotypic traits across both study arms except for autistic features and moderate ID, the latter being associated with a slightly higher, though non-significant, diagnostic yield (31/83, 37.3% vs. 198/567, 34.9% overall) and being more common in the SoC arm (50/283, 17.7% vs. 33/284, 11.6% in the GS arm). As a result, adjusting for this factor would have marginally increased the yield of GS.

The study protocol was approved by the Ethics Committee of UZ Leuven (S64603) after consultation with the ethics committees of all participating sites and was performed in accordance with the Declaration of Helsinki. Given the absence of a therapeutic intervention, the study was not prospectively registered in an external trial registry. However, as participants were randomized to the diagnostic strategies, the study qualifies as an interventional trial following current recommendations from the International Committee of Medical Journal Editors. The trial was therefore retrospectively registered on ClinicalTrials.gov (NCT07051213). The study was performed in accordance with the Declaration of Helsinki. Written informed consent was obtained from all participants or their legal representatives.

Consent for publication was obtained from all participants or their legal representatives.

Illumina Inc. (San Diego, CA, USA) provided the consumables for genome sequencing of the samples used for the ring trial and the 284 probands and their parent(s) sequenced in the genome sequencing arm, as well as access to the Emedgene platform.The authors declare that they have no other competing interests.

Information regarding the phenotype and molecular diagnosis of study participants is shared in the additional files. The dataset generated and analyzed during the current study is available through the EGA repository (EGAC00001003120), https://ega-archive.org/datasets/EGAD50000002905, containing 1166 samples including the samples from centers 1, 2, 4, 5, 8, and part of the GS samples from center 7.