In subject 1 of family 1, we identified a de novo heterozygous missense SNV in the DEAF1 transcription factor gene (MIM 602635) as shown in Table 1. De novo heterozygous mutations in DEAF1 cause intellectual disability, speech impairment, and autistic behaviors (25). Multiple de novo heterozygous missense mutations have been identified in DEAF1 causing Vulto-van Silfhout-de Vries syndrome characterized by ID with severe speech defects and behavioral problems (25–28). Additionally, homozygous recessive mutations in DEAF1 have been identified in subjects with autism, hypotonia, language delay, ID, basal ganglia dysfunction, and epilepsy (29, 30). The variant c.674G > T (NM_021008.4) in exon 5 predicted a p.G225V substitution which was deemed deleterious, because this variant was not reported in gnomAD and highly disruptive by SIFT and Polyphen 2, respectively. The variant had a CADD score of 26.5 and was classified as pathogenic based on ACMG classification (PS2, PS4, PM2, PP2, PP3) and as possibly pathogenic by M-CAP (Table 1) (31). DEAF1 exhibits high levels of expression in the brain and central nervous system, implying its involvement in early neurodevelopmental stages (25, 30).

Within the same subject, we found an additional variant in EPHA5 (EPH receptor A5, MIM 600004), a gene that has not been previously linked to neurodevelopmental disorders. Eph receptors including EPHA5, have been found to modulate the strength of existing synapses in the adult brain (32). The identified variant c.680G > T (NM_182472.5) was homozygous and caused p.R227L change by disrupting exon 3. It was predicted to be likely benign by M-CAP with CADD score 28.5 and classified as likely pathogenic by ACMG due to PS4, PM2, PM3, and PP3 (Table 2).

In family 2, we found two different de novo variants, in chloride voltage-gated channel 3, CLCN3 (MIM 600580) and the ligand dependent nuclear receptor interacting factor 1, LRIF1 (MIM 615354) genes. De novo heterozygous missense variants in CLCN3 were identified in 9 unrelated patients who had neurodevelopmental disorders with hypotonia and brain abnormalities (33). CLCN3 variant c.1952A > C (NM_001243374.2) at 4q33 in exon 10 in subject 2 resulted in p.D651A. This was a deleterious variant with CADD score 22.6 and classified as pathogenic by M-CAP and with ACMG classification of PS2, PS4, PM2, PP2, and PP3 (Table 1).

Another candidate gene LRIF1, located at 1p13.3, has been suggested as a potential neurodevelopmental candidate gene. Its physical interaction with four established neurodevelopmental genes (PQBP1, CHD3, CHD4, and ADHC1) further enhances its appeal as a candidate gene for ASD in this proband (34). The variant in LRIF1 in subject 2 had c.1240A > G (NM_018372.4) change, resulting in p.I414V in exon 2. The CADD score of this variant was 20.6 and based on ACMG classification-PS2, PS4, PM2, PP3, it was denoted as pathogenic, and likely benign by M-CAP prediction (Table 2). A homozygous frameshift mutation in LRIF1 was reported in a man with digenic facioscapulohumeral muscular dystrophy-3, who was born to consanguineous Japanese parents (35).

After considering factors such as the CADD score, M-CAP, and ACMG interpretation, CLCN3 takes precedence over LRIF1 in terms of prioritization.

Within family 3, we identified two de novo variants in TRPC4 (Transient receptor potential cation channel subfamily C member 4, MIM 603651) and RGS22 (Regulator of G protein signaling 22, MIM 615650). Notably, neither of these genes has been previously reported to be associated with neurodevelopmental disorders. TRPC4 located at 13q13.3 harbored a splicing variant c.379-2A > T (NM_016179.4), with a high CADD score of 34, indicating its deleterious nature. All splice sites at the exon/intron boundaries consistently features dinucleotides. GT at the 5′ splice donor site and AG at the 3′ splice acceptor site. Point mutations within these consensus sequences can disrupt splicing event, leading to the generation of an aberrant transcript from the affected gene (36). According to ACMG classification, this variant was classified as pathogenic, based on criteria PVS1, PS2, PS4, PM2, and PP3 (Table 1). Second candidate gene, RGS22 at 8q22.2, on the other hand, is a member of a protein family involved in the regulation of G protein signaling pathway (37). Overexpression of RGS22 has been shown to reduce cell migration and invasive potential in a metastatic esophageal cancer cell line (38). The RGS22 variant in family 3 was c.2386G > A (NM_001286693.2) in exon 17, resulting in p.A796T. This variant was also deemed deleterious with a CADD score of 23.4 and classified as pathogenic according to ACMG classification criteria PS2, PS4, PM2, and PP3 (Table 2). TRPC4 emerged as the prime candidate gene over RGS22, given its deleterious splicing variant and high CADD score of 34.

In family 4, we detected de novo variants in two lysine methyltransferases, KMT2C (Lysine methyltransferase 2C, MIM 606833) at 7q36.1 and KMT5A (Lysine methyltransferase 5A, MIM 607240) at 12q24. While KMT5A represents a novel gene, KMT2C, also known as MLL3, is a well-established gene linked to intellectual disability and autism (39). KMT2C de novo heterozygous mutations have been associated with severe intellectual disability in 3 unrelated patients (40). Knockdown of KMT2C in the drosophila ortholog “trr” impaired short-term memory in the mushroom body of the fly brain (41). The KMT2C variant was c.2573G > T (NM_170606.3) in exon 15, resulting in p.W858L. The variant was predicted to be deleterious and probably damaging due to a high CADD score of 27 and classified as pathogenic according to ACMG criteria PS2, PS4, PM2, PP2, and PP3 (Table 1). KMT5A (also known as SETD8) plays a role in transcriptional regulation, heterochromatin formation, genomic stability, cell cycle progression, and development through its methyltransferase activity (42). The KMT5A variant c.31C > T (NM_020382.7) in exon 2 leading to p.R11C was deemed deleterious with a CADD score of 20.5, and classified as pathogenic according to ACMG criteria PS2, PS4, PM2, and PP3 (Table 2). KMT2C was selected as the primary candidate gene due to its significantly higher CADD score of 27 and its established association as a NDD gene.

Within family 5, we identified a sole candidate variant in SCFD2 (Sec1 family domain containing 2) at 4q12. The homozygous missense variant c.1148C > T (NM_152540.4), p.P383L in the exon 4 was deemed deleterious and probably damaging with a high CADD score of 27 and classified as likely pathogenic based on ACMG classification criteria PS4, PM2, PM3, and PP3 (Table 1). p.P383L variant is located in sec family domain (aa 377-665) and therefore could cause disturbances in the functional properties of SCFD2 protein. Notably, several studies have proposed this gene as a potential susceptibility gene for neurodevelopmental disorders including autism. A genome-wide association study identified SCFD2 as one of the susceptibility gene among 15 highlighted genes on chromosome 4 for schizophrenia, bipolar disorder, and major depressive disorder (43). Additionally, a comprehensive review listing 792 known or proposed ASD susceptibility genes identified SCFD2 as a contributor to ASD (44). Interestingly, another report spotlighted SCFD2 as an autism-implicated gene with high-confidence single nucleotide polymorphisms in its 3’UTR miRNA recognition element, potentially affecting its expression (45).

In family 6, we discovered two de novo and two homozygous recessive variants in four different genes. DHX30 and PDPR were found to have de novo variants, while LRRC18 and PARP14 exhibited autosomal recessive variants. DHX30, known as DExH-box helicase 30 (MIM 616423), encodes an RNA-dependent RNA helicase that utilizes ATP hydrolysis to unwind RNA secondary structures, thereby regulating RNA metabolism and function. Previous research identified six different de novo heterozygous missense mutations in 12 patients with neurodevelopmental disorders along with severe motor impairment and language delay in DHX30 at 3p21.31 (46). The splicing variant c.1813-1G > T (NM_014966.4) in this proband was identified as deleterious, supported by its remarkably high CADD score of 33, suggesting a probable disruption of the protein. It was classified as pathogenic by ACMG classification criteria PVS1, PS2, PS4, PM2, PP2, and PP3 (Table 1).

PARP14, Poly(ADP-ribose) polymerase family member 14 (MIM 610028), plays a role in poly(ADP-ribosylation), that is an immediate DNA damage-dependent post-translational modification of histones and other nuclear proteins, helping in the survival of injured proliferating cells (47). In exon 6 of PARP14 at 3q21.1, we identified the variant c.2906C > G (NM_017554.3). The variant resulted in a deleterious substitution of p.A969G with a CADD score of 22.6. It was classified as possibly pathogenic by M-CAP and likely pathogenic by ACMG classification of PS4, PM2, PM3, PP3 (Table 2).

PDPR (Pyruvate dehydrogenase phosphatase regulatory subunit, MIM 617835) inhibits the activity of the catalytic subunit (PDP1, MIM 605993) in a Mg (2 +)-dependent and Ca (2 +)-stimulated manner (48). In exon 7 of PDPR at 16q22.1, we identified the variant c.697G > C (NM_001322118.1), resulting in deleterious p.E233Q with a CADD score of 25.3. According to ACMG classification PS2, PS4, PM2, and PP3, it was classified as pathogenic (Table 2).

LRRC18 (Leucine rich repeat containing 18, MIM 619002) at 10q11.23, known for its involvement in the regulation of spermatogenesis (49), harbored the variant c.416_417delCAinsTG (NM_001006939.4) in exon 2. It led to a deleterious variant p.T139M with a CADD score of 25.4. Based on ACMG classification (PS4, PM2, PM3, PP3), it was classified as likely pathogenic (Table 2).

To date, PDPR, LRRC18, and PARP14 have not been associated with neurodevelopmental disorders. In family 6, DHX30 was selected as the primary candidate gene over PARP14, PDPR, and LRRC18 due to its significantly higher CADD score of 33, its classification as pathogenic according to ACMG criteria, and its established association as a NDD gene.

In summary, our study has revealed the involvement of four known genes (DEAF1, CLCN3, KMT2C, and DHX30) as well as two novel genes (TRPC4 and SCFD2) in autism spectrum disorder across six families from Qatar.

The effects of missense variants in ten candidate genes were analyzed at an atomic level using molecular modeling. Homology protein models were employed to examine three-dimensional structural models of various proteins. Crystal structures were utilized for EPHA5 and PARP14, while Alphafold predictions from the website https://alphafold.ebi.ac.uk were used for DEAF1, CLCN3, SCFD2, LRRC18, LRIF1, RGS22, KMT5A, and PDPR. For the DEAF1 variant G225V, the structural modeling based on Alphafold prediction revealed a severe steric hindrance between the side chain of V225 and the main chain carbonyl atoms of Arg218 and Ala283. As a result, the G225V substitution is expected to impair proper protein folding and induce structural instability within the SAND domain (residues 193-273) of DEAF1 (Figure 3, top).

In the crystal structure of EPHA5, Arg227 mediated intramolecular electrostatic interactions with Glu60 and Glu129. Additionally, it formed a hydrogen bond with Asn62. The R227L substitution in EPHA5 induced alterations in the surface shape and charge of the ligand-binding domain (LBD; residues 60–238), which might affect the association of EPHA5 with its putative binding partners. However, it is important to note that Arg227 is not located in the ligand-binding surface of the EPHA5 LBD, indicating that this variant does not appear to interfere with its ligand binding (Figure 3).

In CLCN3 (D651A) modeling, no severe alterations at a molecular level were observed by the D651A substitution in CLCN3. However, since Asp651 is located at the cytosolic C-terminal cystathionine β-synthase (CBS) domain, there is a possibility that the D651A substitution could affect the C-terminal domain-mediated protein recognition of CLCN3 (Figure 3).

Regarding SCFD2 (P383L) modeling, severe steric hindrance was observed between the side chain of L383 and the main chain carbonyl of Ala414, as well as the main chain amide of Lys418 in the Alphafold-based structural model. Consequently, the P383L substitution is expected to impair the proper protein folding, leading to structural instability of the SCFD2 (Figure 3). Analysis of crystal structure of PARP14 indicated that the A969G substitution led to the loss of A969-mediated intramolecular hydrophobic interaction with L930 and V934, which should negatively affect the proper protein folding (Figure 3). In the Alphafold-based structural model of LRRC18 (T139M), steric hindrance was also evident between M139 and V115/N118. Therefore, the T139M substitution is likely to disrupt the proper protein folding of the leucine-rich repeats and cause structural instability of LRRC18 (Figure 3). While protein modeling of the PARP14 with A969G variant and the LRRC18 with T139M variant indicated protein misfolding in family 6, their impact is considered minor compared to the potential protein truncation resulting from the splice variant of DHX30.

Structural modeling of mutations in the disordered loop region is challenging due to their high flexibility. This limitation applies to LRIF1 (I414V), KMT2C (W858X), and KMT5A (R11C). It is worth mentioning that isoleucine and valine share similar structural and biochemical characteristics, and therefore the I414V substitution in LRIF1 is not expected to induce significant molecular-level alteration (Supplementary Figure 1). Furthermore, no significant molecular-level alterations were observed for the A796T substitution in RGS22 and the E233Q substitution in PDPR. It is difficult to interpret their biological consequences at the molecular level. Protein modeling for TRPC4 in family 3 and DHX30 in family 6 was not performed due to potential truncation resulting from splice variants.

In the first family, we identified two potential missense candidate variants, one in DEAF1 as a de novo and another in EPHA5 as a homozygous variant. DEAF1 interacts with several genes at the protein level, including FBXO11, DDX58, AIMP2, GSK3B, and HRAS. They are associated with NDDs including autism, ID, DD, etc. Multiple de novo variants in FBXO11 (F-Box only protein 11) have been found to disrupt its protein expression and localization, leading to neurodevelopmental disorder (NDD) through haploinsufficiency (50). AIMP2 has been identified as a candidate gene for late-onset Parkinson’s disease, expressed in the mesencephalon, and associated with dysregulation in PD associated pathways (51). Another interacting protein, HRAS-germline mutations were involved in Costello syndrome, a multiple congenital anomaly and intellectual disability syndrome (52). Since family 1 is not consanguineous, we chose to give priority to the well-known autosomal dominant neurodevelopmental gene DEAF1 over EPHA5, which has a homozygous variant. This was supported by the results of protein modeling, which revealed a significant steric hindrance caused by the DEAF1 variant, whereas the EPHA5 variant did not seem to affect its ligand binding.

In the second family, CLCN3 emerged as the prime candidate gene over LRIF1. CLCN3 is a member of the CLC family of Cl- channels and Cl-/H + exchangers located in endosomal/lysosomal compartments. Loss of function variants in CLCN3 have been associated with a variety of NDDs such as GDD, ID and dysmorphic features (33). LARP7, an interacting protein of CLCN3 at the protein level, has been reported as a causative of Alazami syndrome characterized by failure to thrive, short stature, and DD symptoms (53). Another interacting protein of CLCN3, SLC35A2 has been linked with a rare type of congenital disorder of glycosylation named as SLC35A2-congenital disorders of glycosylation described by neurological impairments with or without skeletal abnormalities (54). Brain somatic SLC35A2 mutations are also linked to neocortical epilepsy (55). After evaluating factors such as the CADD scores, M-CAP, and ACMG interpretation, the prioritization favors CLCN3 over LRIF1. This choice is supported by protein modeling results, which indicate that the CLCN3 variant has the potential to impact protein recognition mediated by the C-terminal domain, while the LRIF variant is not anticipated to cause substantial molecular-level alterations.

In family 3, we identified TRPC4 as a prime candidate gene having a de novo splice variant. TRPC4, transient receptor potential cation channel, subfamily C, member 4, forms homodimers and heterodimers with other family members such as TRPC1 and 5 (56, 57). Multiple sporadic variants in TRPC family members associated with NDDs have been reported (58–60). ITPR1 (inositol 1,4,5-trisphosphate receptor, type 1), an interacting partner of TRPC4, is causative for Gillespie syndrome characterized by bilateral iris hypoplasia, congenital hypotonia, non-progressive ataxia, progressive cerebellar atrophy, and intellectual disability as well as autism and schizophrenia (61, 62). TRPC4 takes precedence as the prime candidate gene over RGS22 due to its deleterious splicing variant and a high CADD score of 34. This decision was supported by protein modeling results, which indicated no notable molecular-level alterations with the RGS22 variant.

In family 4, we found two de novo missense variants in two lysine methyltransferases, KMT2C and KMT5A, respectively. The well-known neurodevelopmental gene KMT2C interacts with multiple proteins, WDR5, CUL3, HCFC1, and SETD1A that are involved in various NDDs. WDR5 controls neuronal migration and dendrite polarity of pyramidal neurons by affecting reelin signaling (63). CUL3 has been associated with different NDDs including ASD, DD, and ID (64, 65). HCFC1 mutations cause intellectual disability by affecting normal brain functioning (66). Dominant de novo LoF variants in SETD1A, a H3K4 methyltransferase, affect normal cognitive functioning via disbalancing neuronal processes and cause NDDs (67). KMT2C was prioritized as the primary candidate gene over KMT5A based on its notably high CADD score 27 and its established association as a neurodevelopmental gene.

In family 5, our analysis identified SCFD2 as the sole candidate gene. According to Biogrid interaction analysis, SCFD2 interacts with multiple protein partners involved in a variety of NDDs. Notably, PCDHGC4, one such partner, is a neurodevelopmental gene associated with progressive microcephaly, seizures, and joint anomalies, (68). Furthermore, bi-allelic mutations in HS2ST1, another SCFD2 interactor (69), have been associated with developmental delay, intellectual disability, corpus callosum agenesis, facial dysmorphism, and skeletal and renal anomalies (70).

In family 6, we identified one known NDD gene DHX30 with a de novo splice variant as the prime candidate gene. It had 0% gnomAD frequency and a CADD score of 33. CHD3, an interacting partner of DHX30, is known to cause a neurodevelopmental syndrome with macrocephaly, impaired speech, and language delay (71). Another interacting partner CHD4 is involved in a multisystemic neurodevelopmental disorder (72). Furthermore, DHX30 forms an interaction with CUL4B, which is associated with X-linked intellectual disability syndrome (73).

To check the importance of the mutant amino acid (P383) in SCFD2, its conservation across eight vertebrate species was checked by multiple alignment. The SCFD2 protein sequence alignment revealed that P383 found mutated in this study was fully conserved among eight vertebrate species namely, human, chimpanzee, macaque, wolf, cattle, mouse, rat and bird. This conservation in protein sequence across eight vertebrate species emphasized the detrimental consequences of this substitution accentuating its vital involvement in shared biological processes across these species (Figure 4).

The subject 4 had a FMR1 full mutation, being positive for Fragile X-syndrome (FXS) testing (> 200 CGG repeats in the 5′-untranslated region of FMR-1 gene) along with a variant in KMT2C, suggesting a putative digenic model (74). There has been a report of concomitant diagnosis of FXS with three other genetic conditions such as Duchenne muscular dystrophy (DMD), PPP2R5D-related and MYT1L-related intellectual disability (74). To date, no prior report has documented the association between an X-linked recessive FMR1 variant and an autosomal dominant KMT2C variant in autism and intellectual disability. With the rising number of cases involving dual genetic diagnoses identified through exome and genome sequencing, it is crucial to recognize the possibility for digenic pattern in the molecular diagnosis of seemingly monogenic disorders.

ES, a cost-effective alternative to GS, detects variants in the coding regions of the genome, known as the exome, which harbor up to 85% disease-causing mutations. The workflow of ES involves isolating genomic DNA from blood samples obtained from affected individuals and their parents as controls. The isolated DNA is then used to prepare exome-enriched libraries using AgilentSureSelectXT kit. The libraries, after enrichment, are pooled into 16 different pools and clustered on a cBot platform. Subsequently, the samples are sequenced on an Illumina HiSeq 4000 using paired end reads (2 bp × 150 bp) at a coverage of 30X per sample. The obtained next generation sequencing data is aligned to the human genome reference sequence, (hg19) and variant are called using the CLC genomics software from Qiagen. Finally, clinically relevant variants are annotated using QC-IT from Qiagen. In this study, we analyzed the ES data obtained from six syndromic autism trios (affected child and both non-affected parents) from Qatar.

Following the manufacturer’s recommendations, the QIAamp DNA blood kit (Qiagen, Hilden, Germany) was used to isolate DNA from peripheral blood samples collected from each participant. Fragile X screening was performed on each patient, and patients with FMR1 gene change were excluded from all future genetic analyses and phenotypic cluster analysis. This exclusion was necessary because Fragile-X syndrome is a well-known cause of syndromic autism and intellectual disability and has a distinct molecular pathomechanism.

We employed a combination of several techniques to assess the pathogenicity of variants identified in the VCF files. The initial list of candidate variants was obtained by implementing the following filtering steps: variants with a frequency below 0.1% in ExAC/gnomAD v2.11/1000g2015 were included; only exonic, splicing, non-synonymous, and stopgain variants with a sufficient coverage were considered; variants predicted to be benign or tolerated in PolyPhen and SIFT were discarded, while those with CADD score ≥20 were retained; variants present in more than 1% of in-house controls (i.e., individuals without ASD or any neuropsychiatric disease) were excluded. Finally, M-CAP analysis also indicated likely pathogenicity of the identified variants (31). Furthermore, we performed additional filtering based on annotation data from public databases such as the Mouse Genome Database (MGD), OMIM, PubMed, and ClinVar. The final section of variants was guided by ACMG guidelines, absence in control datasets (gnomAD and QGP), relevance of interacting proteins in neurodevelopmental disorders and autism based on literature, presence of sporadic variants in the Human Gene Mutation Database (HGMD), expression patterns in human tissues, and phenotypes observed in available knock-out/deletion model organisms. Protein modeling is employed to assess the impact of missense variants on the structure and functionality of our candidate genes identified. By following this rigorous procedure, we obtained a final list of six candidate variants, including five de novo and one recessively inherited.

Molecular modeling of potential candidate genetic variants was conducted using the three-dimensional protein structures obtained from X-ray crystallography [(PDB codes: 4ET7 for EPHA5 and 4D86 for PARP14) or constructed by Alphafold-based prediction]. For the Alphafold protein structures, the corresponding entries in the Alphafold protein structure database² were used: O75398 for DEAF1, P51790 for CLCN3, Q5T3J3 for LRIF1, Q8NE09 for RGS22, Q2YDW7 for KMT5A, Q8WU76 for SCFD2, A8MT40 for PDPR, and Q66HD6 for LRRC18. However, KMT2C was not registered in the database, so the location of its Trp858 in the disordered loop region was predicted based on information from the Uniprot database.³ To introduce the variants into the protein structures and generate visual representations, we utilized the PyMOL program.⁴