Children who visited Shanghai First Maternity and Infant Hospital with suspected motor coordination disorder were assessed using the Movement Assessment Battery for Children-Second Edition (MABC-2). In combination with the DSM-5 (Diagnostic and statistical manual, 2013) diagnostic criteria, the participants were assigned to DCD group according to the diagnostic results. Age-matched healthy children were selected as controls. Peripheral blood samples were collected from all enrolled children for genome-wide DNA methylation analysis. The other 41 participants were from a previous cohort study conducted at the Shanghai First Maternity and Infant Hospital. Motor development at 1 year of age in 41 participants was assessed using Bayley Scales of Infant and Toddler Development, Third Edition (BSID-III), and peripheral blood was collected for MethylTarget™ sequencing. The study was approved by the Ethic Committee of Shanghai First Maternity and Infant Hospital (KS1630). All information acquired was kept confidential and was only accessible by the researchers.

In this study, the diagnosis of DCD was based on both the MABC-2 scores and the criteria outlined in DSM-5. According to the DSM-5 (Diagnostic and Statistical Manual, 2013), DCD should be diagnosed based on the following criteria: 1) acquisition and execution of coordinated motor skills are below the expected level for age, given the opportunity for skill learning; 2) motor skill difficulties significantly interfere with activities of daily living and impact academic/school productivity, prevocational and vocational activities, leisure, and play; 3) onset is in the early developmental period; and 4) motor skill difficulties are not better explained by intellectual or visual impairment or other neurological conditions that affect movement. Additionally, children with co-occurring ASD, ADHD, or learning disorders (LD) were excluded from the study.

The MABC-2 is a widely used diagnostic tool for assessing DCD (Henderson et al., 2007). The test consists of eight items and the recorded raw data are standardized and converted into a standard score of 1–19, reflecting three motor subtests (manual dexterity, aiming and catching, and balance), which represent fine motor skills, gross motor skills, and balance abilities, respectively. Finally, the standardized scores from each subscale are summed to calculate the total MABC-2 score.

The BSID-III is a widely utilized assessment tool for evaluating neurodevelopment in children up to 42 months of age. In this study, we included scores from motor (gross and fine motor) domain. The motor subscale may be useful for describing and assessing motor function, especially for general developmental assessment to identify early motor dysfunction. Frijters et al. (Frijters et al., 2010) have demonstrated that it had a good correlation with MABC-2 results in children aged 36–48 months. As described previously (Hua et al., 2019), the BSID-III was used to assess motor development at 1 year of age in 41 participants. These assessments were conducted by trained professionals. The motor domain composite scores were standardized according to the standard usage of the tool.

Genome-wide DNA methylation analysis was conducted using the Infinium Human Methylation 850K BeadChip (Illumina). Genomic DNA was bisulfite-treated using the EZ DNA Methylation Kit (Zymo Research) according to the manufacturer’s protocol. The treated DNA samples were then hybridized to the BeadChip following the Illumina Infinium HD Methylation Protocol. Raw intensity data (IDAT files) generated from the BeadChip were processed using the ChAMP (Tian et al., 2017) package (version 2.14.0) in R, with the human genome build 19 (hg19) as a reference genome for annotation. DNA methylation levels were represented as β values, ranging from 0 (completely unmethylated) to 1 (fully methylated) for each CpG site. Probes with a detection P value >0.01 and located on the X and Y chromosomes were excluded. Additionally, SNP-related probes and multi-hit probes were removed. To correct for Infinium type I and type II probe bias, the BMIQ (Beta MIxture Quantile dilation) algorithm was applied. The final data set for analysis comprised methylation data from 723838 probes. To assess the presence of variation and potential batch effects in the methylation data, we performed singular value decomposition (SVD) analysis and generated a scree plot to visualize the proportion of variance explained by each principal component, thereby identifying potential sources of variation. Additionally, a quantile-quantile (QQ) plot of observed versus expected–log10(p) values was generated to evaluate the inflation of test statistics and detect any systematic bias. To control for cell type heterogeneity in blood samples, we applied a reference-based cell composition correction using the champ.refbase.fix function in the ChAMP package in R with default parameters, which estimates and adjusts for the proportions of major blood cell types (CD8⁺ T cells, CD4⁺ T cells, natural killer cells, B cells, monocytes, and granulocytes) based on the method by Houseman et al. (Houseman et al., 2012). The estimated proportions of these cell types for each group are summarized in Supplementary Table S1.

MethylTarget™ sequencing (Genesky Biotechnologies Inc., Shanghai, China) was used to assess CpG site methylation levels. Following quality control of genomic DNA, target probes primers were designed and single-site PCR conditions were optimized. The primer sequences for target probes are shown in Supplementary Table S2. Optimized primers were combined into a multiplex PCR panel. After bisulfite conversion, the multiplex PCR panel was used for amplification, ensuring balanced target site products. Indexed primers were then used for PCR to introduce Illumina-compatible tags and sequencing libraries were constructed. Finally, after fragment size verification using Agilent 2100 bioanalyzer, sequencing was performed on the Illumina Hiseq platform (2 × 150 bp) to generate FastQ data.

Statistical analysis was performed using R software (version 4.4.3) and SPSS 26.0 (IBM, Armonk, NY, United States). DMPs between groups were identified using the champ.DMP function in the ChAMP package (version 2.14.0). Given the relatively limited sample size, we applied selection criteria based on both unadjusted statistical significance and effect size. Specifically, DMPs were defined as CpG sites with an unadjusted P value <0.05 and an absolute methylation difference (Δβ) > 0.1. This approach has been successfully applied in previous studies to identify DMPs in small samples (Maltby et al., 2020; Imran et al., 2021; Yao et al., 2021; Maltby et al., 2017). DMRs were identified using the champ.DMR function implemented in the ChAMP R package. We applied two algorithms within this function: Bumphunter and ProbeLasso. For Bumphunter, the minimum number of CpGs per DMR was set to ≥2, and the maximum DMR length was restricted to less than 300 bp. For ProbeLasso, default settings were used as provided by the ChAMP package. The DMRs detected by these two algorithms were compared to evaluate the consistency and robustness of the results across different computational methods.

To assess the genomic distribution of DMPs, we compared the numbers of total, hypermethylated, and hypomethylated DMPs in gene-related and CpG island-related regions with the background distribution of all probes on the Illumina 850K array using Fisher’s exact test. Multiple testing correction was performed using the Benjamini–Hochberg (BH) method and an adjusted P value <0.05 was considered statistically significant. For the chromosomal distribution analysis, a chi-squared goodness-of-fit test was used to compare the observed number of DMPs per chromosome with the expected distribution based on the total number of probes per chromosome on the 850K array.

Functional enrichment analysis, including Gene Ontology (GO) (Ashburner et al., 2000) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2008) pathway analysis, was performed using the Over-Representation Analysis (ORA) method implemented in the R package clusterProfiler (version 3.12.0) (Yu et al., 2012) and the results were visualized using the SRplot platform (www.bioinformatics.com.cn).

In the replication phase, the associations between methylation levels at candidate DMPs and motor outcomes (total motor, gross motor, and fine motor scores) were evaluated using univariate and multiple linear regression models. For the multiple regression models, covariates included maternal age at delivery, child’s gender, and birth weight. Maternal age was grouped into two categories: “<35” and “≥35” years. Effect estimates were presented as standardized β coefficients with 95% confidence intervals (CI). A two-tailed P value <0.05 was considered statistically significant. Sensitivity analysis was conducted by adjusting for the covariates to assess the robustness of the findings.

The diagnosis of DCD was performed according to DSM-5 criteria, and 9 of the tested children had confirmed DCD. One child over 7 years of age was excluded based on age (3–7 years) and the remaining participants were then matched in groups based on age and gender. Of the 16 participants, 8 children with DCD and 8 children without DCD had mean ages of 3.25 ± 0.71 and 3.75 ± 1.39 years, respectively, with no statistically significant difference between the groups (t = 0.907, P = 0.3795). The sex distribution was identical between the DCD and control groups, with each group including 5 boys and 3 girls. The MABC-2 score percentiles were significantly lower in the DCD group (6.30 ± 8.72) compared to the control group (74.38 ± 19.19), with a statistically significant difference (t = 9.134, P < 0.0001). Peripheral blood was collected from enrolled children for methylation data analysis (Table 1).

Genome-wide DNA methylation profiling of 16 children was detected using 850K BeadChip. To assess potential technical batch effects, SVD analysis was conducted using variables including slide, array, plate, well, start_date, and sample_group. As shown in Supplementary Figure S1A, only the sample_group variable was significantly associated with the first principal component (PC1, P < 0.05), whereas all technical variables were not significantly associated (P > 0.05). A scree plot (Supplementary Figure S1B) further demonstrated that the first three principal components explained 41.82% of the total variance. These findings suggested that the primary source of variation was biological rather than technical. Additionally, the QQ plot has been presented in Supplementary Figure S1C. The observed genomic inflation factor (λ = 1.168) falls within acceptable limits, indicating no substantial inflation or systematic technical bias in the test statistics.

Based on these results, we subsequently performed differential methylation analysis between groups. As shown in Figure 1A, a total of 416 DMPs, including 151 hypermethylated DMPs (36.30%) and 265 hypomethylated DMPs (63.70%), were identified in the DCD group compared with the control group. 277 of these DMPs had UCSC RefGene Name, mapping to 238 unique genes, while other DMPs were not annotated to specific genes. The RefGene name and the location of each DMP in a gene or chromosome are provided in Supplementary Table S3 and Supplementary Figure S2. We next examined the genomic distribution of these DMPs. As shown in Figure 1B, DMPs were most frequently located in gene bodies (Body, 38.22%) and intergenic regions (IGR, 33.41%). However, Fisher’s exact test indicated that DMPs were significantly enriched in IGR (adjusted P = 0.047) but not in Body regions. In contrast, DMPs were significantly depleted in TSS200 regions (adjusted P = 0.047) (Figure 1B; Supplementary Table S4). Similarly, within CpG island-related regions, DMPs were significantly enriched in opensea regions (adjusted P = 4.36 × 10⁻⁶) and depleted in CpG islands (adjusted P = 1.04 × 10⁻⁶) (Figure 1C; Supplementary Table S4). Further stratified analysis showed that hypomethylated DMPs were also enriched in opensea regions (adjusted P = 1.56 × 10⁻⁴) and depleted in CpG islands (adjusted P = 4.12 × 10⁻⁶), while no significant enrichment or depletion was observed for hypermethylated DMPs in any genomic or CpG island-related region (Figures 1D,E; Supplementary Tables S5, S6). The distribution of DMPs varied among chromosomes, with the highest number observed on chromosome 1, chromosome 6, and chromosome 2 (Figure 2A; Supplementary Table S7). However, a chi-squared goodness-of-fit test showed that this distribution did not significantly differ from that of all detected probes (χ² = 26.13, df = 21, P = 0.2016), indicating no significant chromosomal enrichment. Additionally, a heatmap showing the methylation levels of all DMPs is presented in Figure 2B. The top 10 DMPs ranked by P value are shown in Figure 2C, and those ranked by absolute methylation difference are shown in Supplementary Figure S3.

A total of 48 DMRs were identified between the DCD and control groups using the Bumphunter algorithm, mapping to 48 unique genes (Supplementary Table S8). The distribution and methylation level of these DMRs were visualized by a circos plot and heatmap (Figures 3A,B).Among these, 12 genes overlapped with those annotated by DMPs. A hypergeometric test showed that the overlap is significantly more than expected by chance (P = 2.91 × 10⁻¹³). To assess the robustness of DMR detection, we further applied the ProbeLasso algorithm implemented in the ChAMP package with default parameters. This analysis identified 22 DMRs mapping to 24 genes (Supplementary Figure S4A; Supplementary Table S9). Among these, three genes (ABAT, NXN, and MIR365A) were consistently identified by both ProbeLasso and Bumphunter (Supplementary Figure S4B), suggesting that these genes represent robust DMR signals across different computational strategies. This cross-validation supports the reliability of the observed regional methylation changes. For subsequent analyses, we used the DMRs identified by the Bumphunter algorithm.

To explore the functional annotation of DMPs and DMRs, we performed GO term and KEGG pathway enrichment analyses on genes associated with DMPs and DMRs using the ORA method implemented in the R package clusterProfiler. The DMPs-associated genes were enriched in 3385 GO terms and 204 KEGG pathways, while the DMRs-associated genes were enriched in 1356 GO terms and 71 KEGG pathways. To further classify the enriched pathways, we performed KEGG pathway classification analysis (Figure 4A; Supplementary Figure S5C). We found that DMPs-associated genes were enriched in 15 KEGG pathways associated with the nervous system and neurodegenerative diseases (Figure 4B), whereas DMRs-associated genes were enriched in 12 such pathways (Figure 4C). Additionally, the top five enriched GO terms in each category (biological process, cellular component, and molecular function) for DMPs- and DMRs-associated genes are presented in Supplementary Figures S5A, S5B.

To investigate the association between DNA methylation and motor performance in children, we analyzed the correlation between the methylation levels of DMPs and MABC-2 scores. A total of 251 DMPs were significantly associated with motor performance. For further validation, we prioritized DMPs based on their statistical significance, functional relevance, and involvement in biological pathways. We identified 42 key DMPs mapped to genes involved in neurodevelopment and motor function (Table 2). For example, LIN28 and ATXN7 have been implicated in neurogenesis and brain development (Wang and Li, 2024; Niewiadomska-Cimicka and Trottier, 2019). IGF1R plays a crucial role in synaptic plasticity and complex cognitive functions (Cardoso et al., 2021). NCK1 is essential for neuronal connectivity and signaling (Fawcett et al., 2007) and TNFAIP3 is involved in regulating microglia activation and neuroinflammation (Voet et al., 2018). Additionally, the COL5A1 gene encodes the α1 chain of type V collagen, which is crucial for musculoskeletal development (Chandrasekaran et al., 2024).

Next, we validated the associations between the methylation levels of key DMPs and motor performance in a larger population using MethylTarget™ sequencing. Table 3 shows the characteristics of the participants in the replication phase. As shown in Figure 5 and Supplementary Material 2, after adjusting for potential confounders, methylation levels at cg18187326 (FAM45A, chr10:120873428) and cg11968956 (FAM184A, chr6:119441387) remained significantly associated with total motor performance (cg18187326, adjusted β = 0.382, P = 0.018; cg11968956, adjusted β = 0.334, P = 0.041). For gross motor, higher methylation levels at cg18187326 and cg11968956 were also positively correlated with better gross motor performance (cg18187326, adjusted β = 0.444, P = 0.005; cg11968956, adjusted β = 0.405, P = 0.011). In addition, cg05986449 (GPD2, chr2:157320849) was significantly associated with improved gross motor performance (adjusted β = 0.330, P = 0.042). Notably, cg03597174 (SEZ6 (dist = 26417), PIPOX (dist = 10043), chr17:27359875) showed a negative association with gross motor performance, suggesting that increased methylation level at cg03597174 may be linked to poorer gross motor performance (adjusted β = −0.351, P = 0.038). However, no significant associations were observed between methylation levels at these DMPs and fine motor performance. The genomic locations and other annotations of these four DMPs are summarized in Table 4.

We conducted sensitivity analysis to examine the robustness of the associations between methylation levels at cg18187326, cg11968956, cg03597174, and cg05986449 and motor performance. Different combinations of covariates, including maternal age at delivery, gender, and birth weight, were tested. The results remained consistent with the main analysis, supporting the stability of the observed associations (Table 5).