The National Institute for Biotechnology and Genetic Engineering collected a sample of families who suffer from intellectual disabilities in Pakistan. The probands’ biological parents were consanguineous. Among the probands and their affected relatives, deficiencies in movement, linguistic competence, and adaptive behavior were commonly observed (Table 1). In addition, some patients exhibited co-occurring symptoms such as epilepsy or skeletal abnormalities. Three peripheral blood samples were drawn from participants who, or whose guardians, provided consent. Informed consent for the study was obtained, and sample collection and information records were authorized. All procedures were authorized by informed consent from study participants. The study was conducted according to the principles outlined in the Declaration of Helsinki, and all protocols, consent forms, and procedures for the collection of samples were approved by the local institutional review board, specifically the Research Ethics Committee of the National Institute for Biotechnology and Genetic Engineering in Faisalabad, Pakistan, on 20 June 2016.

WES was performed for two probands and their father. Peripheral blood white cells were used to extract whole genome DNA. The exomes were enriched and indexed using BGI exon 59M V4 kits. Sequencing was performed as 100 bp paired-end reads on BGISEQ-500 systems, which used combinatorial Probe-Anchor Synthesis and improved DNA Nanoballs technology (Fang et al., 2018). The data from this study have been deposited in the CNGB Nucleotide Sequence Archive and are publicly available (CNSA: https://db.cngb.org/cnsa;accession number CNP0002171).

Adapter sequences were removed, and low-quality reads, defined as those containing excessive Ns or low base quality, were filtered out. The resulting high-quality “clean reads” were aligned to the human reference genome hg19 using both BWA and SOAP alignment tools. Variation detection was processed using the pipeline from the Genome Analysis Toolkit.

All variants were then annotated for transcripts and mutation types. Quality control was an integral part of our pipeline, ensuring that only high-quality data was used throughout the analysis process, including read alignment and variant calling. The SNPs were called by GATK HaplotypeCaller using a joint calling strategy, which combines information from multiple samples to improve accuracy.

Just as the filtering parameters commonly used (Han et al., 2020; Zhang et al., 2020), we use the following parameters to filter variants here: 1) The VAF should be limited to less than 0.005 and must be segregated by heredity. 2) The variants should be located in protein-coding regions. 3) The variants that are classified as benign or likely benign according to the ACMG guidelines (Richards et al., 2015; Oza et al., 2018) should be eliminated, but those of uncertain significance and above could be candidates. 4) The loci should be conserved across mammals by multiple sequence alignments. 5) The variants should be co-segregated with other family members according to the Sanger validation. Twelve variants were conserved across species, and Sanger was used to validate them. Finally, the variant DOP1A(NM_001385863.1):c.2561A>G was co-segregated within the family by Sanger validation (Table 1; Supplementary Figure S1). Both parents were heterozygous for the DOP1A(NM_001385863.1):c.2561A>G, and the proband inherited the variant from both parents, resulting in a homozygote. The presence of the same rare variant in both parents might suggest consanguinity.

The variants were classified according to the American College of Medical Genetics and Genomics (ACMG) guidelines (Richards et al., 2015; Oza et al., 2018) and the Standard Operating Procedure of the Clinical Genome Resource (ClinGen) General Sequence Variant Curation Process, Version 2.0. According to the current consensus, variants designated as “Likely pathogenic” or higher are generally regarded as having potential deleterious effects (Richards et al., 2015; Oza et al., 2018).

The reference protein and DNA sequences for the upstream and downstream regions of the variants in diverse mammalian species were obtained from the National Center for Biotechnology Information (NCBI,https://www.ncbi.nlm.nih.gov/). Multiple sequence alignments were performed using the ClustalW online tool (https://tcoffee.crg.eu/apps/tcoffee/). Conservation of loci across mammalian species was ensured by retaining only conserved regions in alignment.

Genomic DNA was extracted from peripheral blood samples and amplified. The amplified DNA was then sequenced using a 3730XL DNA analyzer platform (Thermo Fisher). For validation purposes, samples were collected from the two affected children and their unaffected siblings, as well as from their parents. The primer sequences used in this experiment are listed in Supplementary Table S1.

The gene expression dataset was obtained from the GTEx portal gene tpm dataset (Version 8, https:// www.gtexportal.org/home/). The expression levels of DOP1A, MBP, PLP1, KIFC2, KIF5B, KIF5C, MON2, and MAG were extracted from the GTEx dataset. The expression patterns were visualized using the R package “ggplot2".

In 2020, Chen et al. investigated transcriptional changes in different brain subregions within a 100-mm diameter around amyloid plaques using spatial transcriptomics in an Alzheimer’s disease (AD) mouse model (Chen et al., 2020). They demonstrated that the co-expression network enriched for myelin and oligodendrocyte genes showed early alterations in spatial transcriptomes. Dop1a expression profiles were compared between wild-type and AD mice. Furthermore, the gene expression dataset for human brains during fetal development and after birth was obtained from the Brainspan public datasets (http://www.brainspan.org/). The Pearson correlation coefficient was calculated for each pair of genes, and the expression patterns were visualized using the R package “circlize".

The prediction was processed using the online predictor tool PredictProtein (Qiu et al., 2020; Bernhofer et al., 2021), with the reference sequence Q5TA12_HUMAN selected based on WES filtering results. The variant site was compared with the predictions of secondary structure and solvent accessibility. Protein binding sites were predicted using a 3D structure-independent method.

The precise conformation of the DOP1A protein in three-dimensional space remains elusive, prompting the utilization of a dependable amino acid sequence-based approach for predicting its protein-protein interactions, as described by Romero-Molina (Romero-Molina et al., 2019).

Three individuals were sequenced, resulting in the mapping of 26415 Mb bases to the target regions for individuals IV-1, IV-2, and III-2, respectively. The mean depth of the target region was 286X while the mean coverage of the target region was 99.8%.

The phenotypes of the affected individuals were collected and resembled the phenotypes of dop1a null mutant rats in the research on vacuole formation (VF) rats (Tanaka et al., 2012; 2014). Additionally, Pelizaeus-Merzbacher disease (PMD) patients and PMD-like patients investigated in Alexander Lossos’s research were compared (Lossos et al., 2015). PMD and PMD-like diseases are rare genetic diseases that affect the central nervous system, specifically the myelin sheath and they were reported to be caused by variants in PLP1/MAG. Most of the probands and affected relatives were found deficient in movement, linguistic competence, and adaptive behavior (Table 2). The consistent symptoms of PMD manifested in the affected children, such as ataxia, spasticity, nystagmus, and psychomotor retardation. However, the two children exhibited high phenotypical heterogeneity, as observed in PMD patients. The elder brother exhibited abnormal joints and limbs due to severe involuntary muscle spasms, while the younger brother showed milder symptoms. Similar symptoms were also observed in VF rats, although many human tests could not be directly performed on rats.

Multiple sequence alignment analysis of DNA and protein reveals that the loci surrounding the variant are highly conserved across mammal species (Supplementary Figure S2). That indicates these loci are critical for gene function. The amino acid substitution in this region is intolerant and might impact protein function.

DOP1A has never been associated with any neuronal developmental diseases in humans. However, some studies revealed that the gene could be a potential disease-causation gene of neuronal developmental dysfunctions. As mentioned above, the dop1a null mutation rats inbreds were characterized by body tremors, absent-like seizure, abnormal vacuoles in the red nucleus of the midbrain, the reticular formation in the brain stem, and the white matter of the cerebellum and spinal cord, as well as spongy degeneration in the central nervous system which indicates that the gene is a potential target of neuron disorders (Tanaka et al., 2012; 2014). There has been strong evidence that dop1a null mutations cause abnormal myelinogenesis, and myelinogenesis is universally considered the primary cause of many neurodevelopmental disorders. Therefore, the gene is relevant to phenotype-associated mechanisms and manifestations. All the evidence for variant classification could be applied to the individual variant-gene-disease interpretation (Figure 1).

According to the gnomeAD exomes database, there was a single allele with no observed homozygotes and no occurrences in the gnomAD genomes database, indicating an extremely low population frequency of 0.00000479. The variant was the most prevalent in South Asia, with a population frequency of 0.0000696. The homozygous variant was determined to provide moderate PM2 evidence for an extremely low frequency of recessive hereditary diseases.

The proband’s father and mother each carry one copy of the variant. For recessive disorders, rare homozygous variants in trans could get no more than 0.5 points according to the ClinGen General Sequence Variant Curation Process Standard Operating Procedure to avoid overclassification. Given the extremely low population frequency of 0.00000479 and the presence of two affected siblings in the family, a PM3_supporting evidence with a score of 0.5 was assigned.

In the co-segregation testing, one affected segregation and one unaffected segregation were identified (parents and the first proband were excluded). Additionally, three relatives in the aunt’s family underwent co-segregation testing, but none carried the variant (Figure 2). Notably, a male cousin died in his infancy with no records of ataxia, spasticity, epilepsy, nystagmus, or psychomotor retardation. None of the dop1a null rats showed any symptoms of early death (email confirmed with the author). Therefore, it was inferred that the early death of the male cousin may have been caused by another disease. The LOD ratio was calculated with the following equation for autosomal recessive segregations: Z L o d   s c o r e s = ⁡ log ⁡ 10 1 0.25 a f f e c t e d   s e g r e g a t i o n s × 0.7 5 u n a f f e c t e d   s e g r e g a t i o n s . The Lod score was calculated at 1.102, which provided PP1_supporting evidence (Oza et al., 2018).

There are no pathogenic or likely pathogenic missense variants of DOP1A in the Clinvar database, likely because DOP1A has not been implicated in any Mendelian inheritance diseases. However, the gnomAD database reveals 919 observed missense variants in DOP1A, while the expected number is 1249.5, resulting in a Z score of 3.32. This indicates intolerance to missense changes in DOP1A. Thus, the calculated Z score is 3.32. It indicates that DOP1A is more likely to be intolerant of missense changes, and missense could be a potential pathogenic mechanism. According to the Clingen criteria (Z score >3.09), these findings support the idea that missense variants in DOP1A have a low rate of benign missense variation and that they could be a common mechanism of disease. Additionally, PP2_supporting evidence also supports the variant.

Thus, one PM2_moderate evidence, one PM3_supporting evidence, one PP1_supporting evidence, and one PP2_supporting evidence compose a likely pathogenic variant. The classification supports the high probability that it represents the disease-causing mutation in the affected family.

The variant resulted in a substitution of cysteine for tyrosine at the affected amino acid position, which alters hydrophilicity prediction at this locus. As no crystal or 3D structure of any isoforms of protein DOP1A is available, 3D structure-independent prediction was employed. The prediction indicates the loss of several protein binding sites due to amino acid substitution. The mutant site is moderately conserved and located in a predicted helix region. Although the topological transmembrane structure remains largely unaltered by the variant, the affected locus undergoes a transition from being buried to exposed. Notably, the mutated protein gains an additional N-myristoylation site (GNLRCI) at positions 859–864.

In the study conducted by Lei Lu et al., it was observed that the proteins DOP1A and MON2 interact to form a complex that localizes to various cellular organelles, including the Golgi, endolysosome, and endoplasmic reticulum exit site (Mahajan et al., 2019). The N-terminus of DOP1A was found to interact with kinesin-1 (Figure 3). This could result in the disruption of binding sites within the MEC interaction region, further impacting the function of the DOP1A-MON2 complex.

There are no known protein-protein interactions (PPIs) between PLP/MAG and DOP1A in public PPI databases, and the PPI prediction tool utilized in this study also confirmed the absence of any interactions between PLP/MAG and DOP1A (Figure 4). However, the calculated probabilities of online PPI predictions have identified PPIs between DOP1A-MON2, DOP1A-KLC2, and DOP1A-KIF5B, which were validated by Lei Lu et al.‘s laboratory experiments. The consistency between predictions and laboratory experiments supports the credibility of the method used in our study. In addition, the predicted probability of PPI interaction between DOP1A-MAG was 0.843, whereas the probabilities for MON2-MAG and MON2-PLP interactions were estimated to be 0.536 and 0.887, respectively. These results suggest the high probability that PLP and MAG proteins should be transported by the DOP1A-MON2 complex, and functional variants in DOP1A/MON2 could affect PLP and MAG transportation. The results of the PPI prediction could partially explain the consistency between the phenotypes of the DOP1A mutant family and those of PMD/PMD-like patients caused by PLP variants.

In contrast to the AD mouse model, the normalized DOP1A expression levels of wild-type adult mice were relatively lower (Figure 5). The normalized expression levels were compared between brain subregions of wild-type and Alzheimer’s model mice, and five subregions were identified significant differences by the Wilcoxon test. The expression levels were higher in most of the 5 AD mice brain subregions, while they were slightly lower in HPs_CA1_sp and PTL(median_CNU_WT = 2.982974, median_CNU_AD = 4.96435, p _CNU = 0.0474; median_HPd_CA_slm_WT = 2.982974, median_HPd_CA_slm_AD = 3.730511, p _HPd_CA_slm = 0.01832; median_HPs_CA1_SP_WT = 5.331434, median_HPs_CA1_SP_AD = 4.620066, p _HPs_CA1_SP = 0.004408; median_PTL_WT = 4.951096, median_PTL_AD = 4.800058, p _PTL = 0.03138; median_PERI_WT = 2.982974, median_PERI_AD = 4.812775, p _PERI = 0.03937). It has been reported that the myelination in ApoE4 AD patients had been altered by cholesterol_ anomalous oligodendrocytes (Blanchard et al., 2022) and myelin damage was also observed in the AD dementia development process (Dean et al., 2017). The differentiated expression profiles of dop1a in AD and wild-type mice indicate the potential correlations between dop1a and myelin damage.

To expand upon the relationship between DOP1A and myelination, we conducted an analysis of the expression profiles of DOP1A/MON2 and established myelin protein markers, including PLP1/MBP/MAG, utilizing Brainspan datasets (Supplementary Figure S3). Through co-expression analysis of human brain samples across different developmental stages, distinct patterns were identified. During early pregnancy (8–12pcw), the DOP1A-MON2 complex expressed highly, while the three myelin protein markers expressed relatively less. Our results suggested that PLP1 and MBP exhibit mild negative correlations with the DOP1A-MON2 complex, while MAG showed no correlation. (r_DOP1A_PLP1 = −0.4,p_DOP1A_PLP1 = 0; r_MON2_PLP1 = −0.47,p_MON2_PLP1 = 0; r_MON2_MBP = −0.49,p_MON2_MBP = 0). In late pregnancy, the three proteins positively correlated with MON2 strongly (r_MON2_PLP1 = 0.56,p_MON2_PLP1 = 0.02; r_MON2_MBP = 0.64,p_MON2_MBP = 0; r_MON2_MAG = 0.67,p_MON2_MAG = 0). At critical stages, the correlations of PLP1 with the complex disappeared, while the correlations of MBP reached a severe level (2–3_years old: r_DOP1A_MBP = −0.65,p_DOP1A_MBP = 0; 13–23_years old r_DOP1A_MBP = −0.77,p_DOP1A_MBP = 0) (Figure 6). The co-expressions between DOP1A and the myelin protein markers PLP1/MBP/MAG throughout the lifecycle, indicated that DOP1A, as well as DOP1A-MON2 complex, could play a potential role in the myelinogenesis process. The DOP1A/MON2/KIF5B/KLC2/TPR universally expressed in almost every tissue, with the highest expression level in the brain, while the MBP/PLP1/MAG is dominantly highly expressed in the brain and nerve but has little expression in other tissues. The two different expression profiles not only indicated the physiological association between the DOP1A-MON2 complex and MBP/PLP1/MAG, but also suggested multiple transport mechanisms in myelin proteins.