Epigenetic aging signatures and age prediction in human skeletal muscle

Discovery of aging-associated CpG markers from skeletal muscle

Among the 103 PM muscle samples, 23 male samples aged 18 and 35–78 years were examined using the EPIC array to identify aging-associated CpG markers for Next Generation Sequencing (NGS) and Single Base Extension (SBE)-based model construction (Figure 1A and Supplementary Figure 1). Probes with detection p-value greater than 0.05 were excluded, leaving a total of 91,899 high-quality CpG probes for analysis (Supplementary Table 1). Of these, 105 CpG sites met the threshold of an FDR-adjusted p-value < 0.1. Initially, a two-way hierarchical clustering heatmap of 500 randomly selected CpGs based on linear regression was generated, revealing distinct age-related methylation patterns (Figure 1B). In particular, among CpGs with a regression coefficient (β₁) > 0.005, accelerated sites were 1.5 times more prevalent than decelerated ones (181 vs. 122) (Figure 1C). Overall, 70.7% (64,995 CpGs) showed age-related methylation acceleration, while 29.3% (26,904 CpGs) showed deceleration, which aligns with previous findings [21]. Next, enrichment analysis based on CpG island and genomic position was performed using the Chi-square test. Examination of probe counts revealed that accelerated CpGs were widespread across genomic regions classified by CpG density and transcriptional activity, compared to decelerated sites (Figure 1D). Regarding proportion, 30.8% of accelerated CpGs and 21.7% of decelerated CpGs were identified within transcriptional regulatory sites, such as the transcription start site (TSS), including TSS1500 and TSS200, as well as the 5′ untranslated region (UTR) (Figure 1D).

Figure 1. Identification and characterization of aging-related skeletal muscle-specific CpGs. (A) Schematic overview of experimental workflow. The diagram, created using BioRender, illustrates the workflow starting with the collection of PM muscle tissue, indicated an arrow pointing to the anatomical site clearly demarcated through shaded surrounding areas. DNA was bisulfite-converted and analyzed using the Infinium Methylation EPIC array to identify age-associated CpGs. Cryptic and target CpGs were examined via NGS, while only target CpGs were used for SBE. (B) Two-way hierarchical clustering heatmap of methylation profiles from 23 male PM muscle samples profiled by the EPIC array. Samples are color-coded by age group. The heatmap displays 500 randomly selected CpGs from 91,899 significant age-related CpGs (p < 0.05), with beta values normalized to Z-scores (yellow for higher, blue for lower values). (C) Volcano plot showing regression coefficients (β₁) from age-related linear regression. Each dot represents a CpG with p < 0.05. Yellow and blue dots indicate CpGs with methylation gain and loss with age, respectively. (D) Enrichment of significant CpGs by genomic context. The Venn diagram shows the distribution of methylation accelerated (yellow) and decelerated (blue) CpGs among 91,899 significant CpGs. Bar plots display CpG distributions relative to CpG islands (left bar graph) and chromatin regions (right bar graph) using raw counts (the number of CpG sites) and the Chi-square (χ²) test. Statistical significance is denoted as ^*p < 0.05, ^**p < 0.01 and ^**p < 0.001. (E) Filtering of CpGs through multi-step criteria. The Venn diagram displays the overlap of CpGs filtered by effect size (Δβ ≥ 0.2), model fit (R² > 0.65), and prior muscle-specific clocks (MEAT v.1 and v.2). Twelve and eight final CpGs (underlined) were selected for NGS- and SBE-based modeling. (F) Heatmap of the final 20 CpGs selected for age prediction, grouped by age. Beta values are Z-score normalized; age groups are color-coded.

Among the 91,899 filtered CpGs, 12 CpGs were selected as candidate age-associated markers according to the criteria of a maximum to minimum beta value difference (Δβ) > 0.2, R² > 0.65 from linear regression between age and beta value, and an FDR-adjusted p-value < 0.08 (Figure 1E and Table 1). To improve marker diversity and generalizability, eight additional CpGs overlapping with MEAT v.1 and v.2 were selected from 24,050 sites with Δβ > 0.2, despite R² < 0.65, enhancing applicability across various populations and muscle types (Figure 1E and Table 1). The final set of 20 muscle-specific CpGs, labeled MA_01 to MA_20, represents the Muscle Age numbering (Table 1). Of these, 10 markers displayed age-related methylation acceleration, while the remaining 10 showed deceleration as determined from the EPIC array data (Figure 1F). We compared the 20 markers with probes from GEO datasets to determine if they were previously reported or novel. Methylation data from these markers were compared with GSE50498, which contains VL muscle samples from 24 young and 24 elderly individuals [15], and GSE114763, consisting of VL muscle samples from 8 Europeans [22] (Supplementary Figure 2). Nine markers showed similar methylation patterns, with R² differences below 0.4 between VL muscle data from Caucasian donors and PM muscle data from Korean samples. The remaining markers, however, showed marked discrepancies. These results suggest that while methylation patterns are generally consistent, differences in predictive accuracy may reflect anatomical and population variations. This highlights the need for new, robust age-related CpG markers and a better understanding of their genomic context.

Genomic analysis and functional implications of the 20 aging-associated CpG markers

Among the target 20 CpGs, 10 markers exhibited methylation acceleration. Nine of these were mapped to genes including ZNF549, MNX1-AS1, KIF15, MYCNUT, CDKN2B, FOXS1, IL18BP, RUFY3, and PIPOX (Table 1), while one marker (MA_12) was not associated with any annotated gene. Conversely, the 10 CpGs showing methylation deceleration were linked to nine genes including SLC44A2, CFAP74, TPM3, TWF2, ACTA1, CAND2, SLPI, CLCN1, and DHDPSL. The biological functions and age-associated disease relevance of these differentially methylated genes were evaluated using the Ingenuity Pathway Analysis (IPA) platform to reveal distinct biological pathways influenced by methylation changes with aging.

Functional analysis identified statistically significant differences across 21 canonical pathways (p < 0.05), with the top eight pathways (p < 0.01) predominantly influenced by five key molecules: ACTA1, TPM3, HOGA1, PIPOX, and KIF15 (Figure 2A). These molecules were linked to muscle functions, contraction, and metabolism. Furthermore, 276 diseases associated with transcriptional changes driven by differential methylation (p < 0.05) were identified, with 23 showing stronger significance (p < 0.001), including 13 directly related to myopathy (Figure 2B). Notably, three muscle-specific genes, ACTA1, TPM3, and CLCN1, were associated with hereditary skeletal myopathy and congenital myopathy.

Figure 2. Functional enrichment and transcriptional analysis of aging-associated skeletal muscle-specific CpG loci. (A, B) Ingenuity canonical pathways of aging-related CpGs and classification of diseases and functions. Canonical pathways are ranked by -log(p-value) (right x-axis) and associated molecules shown (A). Disease and function categories are similarly ranked and visualized (B). Pathways directly linked to muscle functions are marked with an asterisk (*). These pathways were identified through IPA program. (C) mRNA Transcriptional Levels of Genes Associated with Aging-Related CpGs. Nine genes showed methylation deceleration at the target CpG sites with aging, while nine others displayed methylation acceleration. mRNA transcription levels of the methylation-decelerated genes are presented in the left graph, and the methylation-accelerated genes are shown in the right graph. All mRNA levels correspond to genes with muscle-specific CpGs associated with aging, measured using RT-qPCR in pectoralis major muscle. Tissue samples from 11 individuals were used and categorized into two age groups. The young group (n = 6) consisted of individuals aged 21, 24, 27, 29, 33, and 35 years, including three males and three females. The old group (n = 5) included individuals aged 63, 65, 68, 74, and 77 years, comprising four males and one female. Genes with target CpGs located in the regulatory regions TSS1500, TSS200, and 5'UTR are marked with a hash (#). Statistical significance was assessed using Student’s t-test, with significant differences indicated as ^*p < 0.05, ^**p < 0.01, and ^***p < 0.001. Data are presented as the mean ± standard error of the mean (SEM), highlighting significant differences across groups. IGR, intergenic region.

To examine how methylation affects transcription, RNA expression of the 18 genes associated with the 20 target CpGs were analyzed using RT-qPCR (Figure 2C). Significant differences between young and old individuals were observed for 14 of these genes. In particular, expression of the CLCN1, IL18BP, and RUFY3 genes exhibited more than two-fold changes, with highly significant statistical differences (p < 0.001). Notably, with advancing aging, mRNA levels of genes exhibiting methylation change in regulatory regions, including the 5′UTR and TSS positions, displayed a positive correlation, except for the SLPI gene (Table 1). This trend is potentially explained by DNA methylation occurring within regulatory regions characterized by low CpG density [23]. The RNA expression findings for these 20 genes are consistent with previous transcriptome and proteome investigations [24–26] (Table 1).

Construction and validation of age prediction model using NGS data

To develop an age prediction model representative of muscle aging, we analyzed NGS data from 20 selected CpG markers. MiSeq sequencing generated approximately 39 million reads that passed the initial quality assessment, with an average of 249,133 read pairs per sample, ranging from 95,134 to 608,471. Following processing with CLC software, the average read count per sample was reduced to 50,892, ranging from 18,043 to 130,532 (Supplementary Figure 3A). Among the 20 target amplicons, 14 demonstrated average coverages between 1,276 and 8,967 reads per CpG, exceeding the 1,000 read threshold for reliable methylation quantification (Supplementary Figure 3B). Although the other six amplicons satisfied the theoretical coverage required for a 90% confidence interval [27], four amplicons with coverage below 700 reads per CpG were excluded from model construction to maintain data reliability.

An overview of the machine learning workflow based on NGS and SBE data is illustrated in Figure 3A. In addition to the 16 primary CpGs from the EPIC array, adjacent cryptic CpGs were detected and considered as supplementary candidate loci for model development, potentially capturing subtle but informative age-associated methylation changes (Figure 1A). Model construction included all 16 amplicons, comprising the 16 target CpGs and 54 additional cryptic CpGs. A total of nine models were constructed using various machine learning approaches. Each model was fitted using a training dataset (n = 71) and evaluated on an independent test set (n = 32), resulting in mean absolute error (MAE) values ranging from 5.537 to 7.017 in the test set (Figure 3B). The detailed formulas of all models are presented in Supplementary Table 2. Of these, the Elastic-net (Ela) model, optimized with an alpha of 0.0178 and an L1 ratio of 0.9852, selected 9 target CpGs and 12 cryptic CpGs across 12 amplicons (Table 2). This model achieved a MAE of 3.184 in the training set and 5.537 in the test set. Although the performance was robust, the discrepancy between training and test MAE implies a risk of overfitting. To address this, we attempted to reduce overfitting by applying stronger regularization and post hoc feature selection, however, these approaches did not improve model performance (Supplementary Table 3). To further evaluate reliability, Leave-One-Out Cross-Validation (LOOCV) was performed on the training set, yielding a MAE of 3.916. Additionally, we performed 5-fold and 10-fold cross-validation, which resulted in MAE values of 3.8629 and 3.8743, respectively. To further assess robustness and determine whether the Ela model captures true age-associated signals rather than random noise, we performed permutation testing, which demonstrated that the model is not simply overfitting to random structure (Supplementary Figure 4). The Elastic Net model was identified as the best-performing approach, exhibiting the lowest MAE in the test set and reduced overfitting compared with other models. Remarkably, the Ela model achieved strong predictive accuracy using only 21 CpGs across 11 amplicons, highlighting the efficiency of a compact marker panel. Collectively, these results establish the Ela model with 21 CpGs as a practical and robust tool for predicting chronological age (Figure 3C and Table 2).

Figure 3. Development and validation of age prediction models using the NGS system. (A) Schematic of the data processing workflow for NGS and SBE-based models, created using BioRender. Data were divided into training and test sets at a 7:3 ratio. Supervised feature selection involved correlation analysis and multicollinearity filtering, followed by model construction using forward selection, stepwise selection, penalized regression, and decision tree models with optimized parameters. Model validation was performed through LOOCV, with performance evaluated by Pearson’s r, R², MAE, and RMSE. Finally, the validated prediction model was applied to test sets. (B) Performance heatmap of nine machine learning models developed from NGS data. Models differ by algorithm (LR, SLR, Ela, Las, Rid, RF, GB, XGB) and CpG set. Metrics (r, R², MAE, RMSE) were calculated for both training and test sets (n = 103). LOOCV was used for training validation. Heatmap color scale reflects relative performance for each column criterion and the top NGS model is indicated in bold. (C) Prediction accuracy of the best-performing NGS model. The best model’s accuracy was evaluated by comparing predicted age with chronological age for 103 samples in both the training (orange) and test (blue) sets (left plot). Residuals between DNA methylation (DNAm) age and chronological age are plotted (right plot). Regression lines are shown with each color, with 95% confidence intervals shaded. Model performance metrics are summarized in the accompanying table. Abbreviations: LOOCV: Leave one out cross validation; LR: Linear regression; SLR: Stepwise Linear Regression; Ela: ElasticNet; Las: Lasso regression; Rid: Ridge regression; XGB: XGBoost; RF: RandomForest; GB: GradientBoosting.

Construction and validation of age prediction model using multiplex SBE data

To ensure practical applicability, we aimed to develop a multiplex SBE system. Initially, we generated methylation data from 68 samples using four multiplexed sets on the SBE platform (Supplementary Figure 5). The association between age and methylation levels at each marker was evaluated using Pearson’s r from linear regression analysis. Among the 20 target CpG markers identified in this study, 3 demonstrated strong correlations (r > 0.7), including MA_05 (r = 0.7488), MA_10 (r = 0.8174), and MA_19 (r = 0.7865), and 12 exhibited moderate correlations (0.5 < r < 0.7) (Supplementary Figure 6). Before constructing predictive models, the methylation data were evaluated for adherence to statistical assumptions and normality. Values for skewness (−0.059) and kurtosis (−0.476) indicated a normal distribution (Supplementary Figure 7A). Normality of residuals was confirmed through a Q-Q plot, and the Durbin-Watson test (1.684) supported the independence of errors without autocorrelation (Supplementary Figure 7B). Through the multicollinearity test, marker MA_14 was excluded due to a variance inflation factor of 11.25, exceeding the threshold of 10 (Supplementary Figure 7C).

To identify the most significant markers, 10 distinct models were developed using the remaining 19 CpGs with various machine learning techniques. All models achieved strong performance, with r > 0.902 and R² > 0.754 in training set, and r > 0.824 and R² > 0.623 in the test set (Figure 4A). LOOCV applied to the training set resulted in MAEs ranging from 3.231 to 6.197 years, while the test set showed MAEs between 5.161 and 9.120 years. Stepwise linear regression (SLR) identified a 7 CpG subset, overlapping with the top-performing 21CpG_Ela model derived from NGS data. The SLR model provided the best balance of predictive accuracy and practical simplicity, achieving an MAE of 5.463 and RMSE of 6.894 in the test set (Figure 4A). Three of these 7 CpGs (MA_13, MA_18, MA_19) coincided with previously identified MEAT markers, while the other four (MA_01, MA_04, MA_08, MA_10) were novel to this study.

Figure 4. Development and validation of age prediction models using the SBE system. (A) Performance heatmap of 10 machine learning models based on the SBE platform. Models differ by algorithm (LR, SLR, Ela, Las, Rid, RF, GB, XGB) and number of CpGs used. The performance metrics such as r, R², MAE, and RMSE values were calculated for both training and test sets. A total of 68 samples were used, with 43 assigned to the training set and 25 to the test set. An additional 35 samples were later included, resulting in a final dataset of 103 samples, consisting of 71 samples for training and 32 for testing, which were used for the final modeling. The heatmap displays relative values using a red-to-blue gradient for each column criterion. The top-performing SBE model is shown in bold. (B) Electropherogram of a unified SBE system using seven CpGs selected from both NGS and SBE models. Capillary electrophoresis was performed using a 3500 genetic analyzer (upper plot). The log₁₀(RFU) values for each peak are shown in a box plot, displaying the mean ± SEM (lower plot). Blue fluorescence represents Guanine (G), representing methylated cytosine, while green fluorescence indicates Adenine (A), representing unmethylated cytosine. (C) Age prediction performance of the best-performing SBE model. The best model’s accuracy was assessed by comparing predicted age with chronological age for 103 individuals in the training and test sets, represented by orange and blue points, respectively (left plot). Residuals between DNAm age and chronological age were plotted (right plot). Regression lines and 95% confidence intervals are shown. Performance metrics for the model are summarized in the accompanying table.

To ensure accuracy and usability, a unified SBE system was established, targeting only the 7 CpGs (Figure 4B). Utilizing this consolidated system, methylation data from 103 samples were collected, and split into a training set (70%, n = 71) and a test set (30%, n = 32). The SLR_7CpG model underwent additional optimization, improving predictive accuracy with r values of 0.951 for both sets. The refined model achieved MAEs of 4.018 years in the training set and 3.797 years in the test set, demonstrating strong robustness for age prediction (Figure 4C). The coefficients and intercept for the final model are presented in Table 2.

Assessment of histological applicability and model validity

To confirm the model’s accuracy and versatility, we expanded the evaluation to include different muscle tissue types. Since tissue-specific cellular composition leads to unique methylation patterns, the model’s performance was tested on cardiac and smooth muscle samples from the heart and uterus, respectively. When the SBE-based skeletal muscle age prediction system (SLR_7CpG) was applied, the heart yielded a high MAE of approximately 24.3 years, indicating poor performance (Figure 5). Similarly, the uterus demonstrated limited accuracy with a MAE of approximately 16.2 years (Figure 5). To assess tissue specificity, EPIC array data were analyzed using linear regression to examine the relationship between chronological age and beta values for the 20 skeletal muscle-specific CpG markers in cardiac and PM muscle tissues. Except for MA_10 and MA_19, the R² values revealed substantial differences between cardiac and skeletal muscle types ranging from 0.261 to 0.683 (Supplementary Figure 8). These findings confirm that the 20 CpG markers and the prediction model are highly specific to skeletal muscle. The diminished accuracy in cardiac and smooth muscle likely arises from differences in cellular composition, such as myocardial cells in the heart and uterine myocytes in the myometrium, which create distinct methylation profiles.

Figure 5. Applicability of the best-performing SBE model to cardiac and smooth muscle tissues. Plots compare DNAm age and chronological age using the best SBE model across tissue types. Left: skeletal (gray) vs. cardiac muscle (red, n = 19); right: skeletal vs. smooth muscle from uterus (purple, n = 19). Skeletal muscle samples (n = 103) are shown in gray on both plots. Solid lines represent the regression relationships for each tissue type, with 95% confidence intervals shaded in corresponding colors. Performance metrics for each tissue are summarized in the accompanying table.

To evaluate the robustness of the assay, we examined the effect of varying DNA input amounts on methylation accuracy and age prediction. Methylation values were calculated from the ratio of methylated to unmethylated sequences generated via PCR amplification. The methylation consistency of the 7 core CpGs included in both SBE and NGS models was assessed using artificially methylated DNA controls (Supplementary Figure 9). The average R2 values were 0.9922 for NGS and 0.9871 for SBE, demonstrating minimal PCR bias. Importantly, when DNA inputs were at least 1.67 ng, the median methylation differences remained under 5% for all markers, except MA_04 in the SBE system (Supplementary Figure 10A). However, MA_13 and MA_18 in the NGS system, as well as MA_19 in the SBE system, showed comparatively greater variability, with standard error of the mean (SEM) values exceeding 0.01 in methylation differences at the 1.67 ng input level. Consistent with previous reports, DNA inputs of at least 10 ng showed reliable performance [28]. Reducing DNA input from 13.33 ng to 6.67 ng increased the MAE by 1.46 years in the SBE model and 2.61 years in the NGS model, both of which were constructed using the same set of 7 CpG markers (Supplementary Figure 10B). Collectively, these results highlight the importance of controlling DNA input to maintain prediction accuracy, with MAEs remaining below 6 years.

Study samples

Ethical approval for this study was obtained from the Institutional Review Board of Seoul National University Hospital (IRB NO. C–1912–053-1087). PM muscle samples were collected via autopsy from 103 deceased human individuals of South Korean descent, aged 18 to 85 years, comprising 80 males and 23 females (Supplementary Figure 1). Cardiac (heart) and smooth muscle (uterus) samples were also collected from the 23 female donors. To minimize the effect of post-mortem interval (PMI) on DNA integrity and methylation status, samples were isolated from non-lesioned tissue areas of uncompromised tissues to ensure integrity with PMIs ranging from 1 to 12 days, averaging 2.4 days. After excision, all tissue samples were promptly stored at −80°C to preserve their molecular stability until further processing.

Genomic DNA extraction

Genomic DNA (gDNA) was extracted from tissue samples using the QIAamp^® DNA Mini Kit (Qiagen, Hilden, Germany). Key adjustments to the protocol included processing 25 mg of tissue with a 4-hour lysis step, optimized to maximize extraction efficiency. The purified gDNA was quantified using the Qubit^™ dsDNA HS Assay Kit on a Qubit^™ Flex Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and subsequently stored at −20°C.

Identification of skeletal muscle-specific aging-related CpG sites

To identify age-associated CpG sites in skeletal muscle, DNA methylation profiles from 23 male PM muscle samples were analyzed, including twenty from the GSE244996 dataset and three newly generated for this study (GSE294234). Genome-wide DNA methylation profiling was conducted using the Illumina Infinium Methylation EPIC BeadChip array version 1 (Illumina, San Diego, CA, USA) in collaboration with Macrogen Inc. Bisulfite conversion was performed on 250 ng of gDNA to facilitate methylation-specific analysis, followed by amplification to ensure sufficient DNA quantity for hybridization. The DNA was enzymatically fragmented, hybridized to sequence-specific probes on the array, and underwent single-base extension incorporating labeled nucleotides. After staining, arrays were scanned on an Illumina scanner. Fluorescent intensity data were extracted and preprocessed in R with the lumi package to perform background correction and dye bias equalization. Probes with detection p-value ≥ 0.05 in more than 25% of samples or missing values (NA) in any sample were excluded. Beta Mixture Quantile (BMIQ) normalization was implemented to standardize methylation values. The resulting Beta and M-values were calculated, and CpG site selection was further refined using statistical criteria, including detection p-value, a maximum beta value difference (Δβ) and R² values from regression between beta values and chronological age. For comparison, VL muscle data from GSE114763 (n = 8) and GSE50498 (n = 47), and cardiac muscle data from GSE244996 (n = 20), were obtained from the Gene Expression Omnibus (GEO).

Genomic analysis of aging-associated genes

To identify enriched signaling and metabolic pathways, as well as to predict the activation or inhibition of upstream regulators and evaluate downstream consequences on diseases, biological functions, and phenotypes, the IPA plugin version 23.0 (Qiagen, Frederick, MD, USA) was utilized. Inferences for expression and pathway analysis were made indirectly using raw p-values and regression coefficients obtained from EPIC array data. Statistical significance was determined by a −log₁₀ (p-value) > 1.3, corresponding to a p-value < 0.05.

mRNA expression analysis by quantitative real-time PCR

Total mRNA was extracted from 10 mg of muscle tissue using QIAzol and the RNeasy Plus Universal Mini Kit (Qiagen) and quantified using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Relative mRNA expression levels were analyzed through Real-Time quantitative PCR (RT-qPCR) on a CFX Connect Real-Time PCR System (Bio-Rad Laboratories, Hercules, CA, USA), utilizing the SensiFAST SYBR Lo-ROX One-Step Kit (Meridian Bioscience, Cincinnati, OH, USA) and custom-designed primer sets (Supplementary Table 5). To better assess the correlation between methylation of the target CpGs and mRNA transcription levels, custom primers were designed to amplify mRNA regions encompassing or adjacent to the 20 CpG sites of interest.

Bisulfite conversion

Bisulfite conversion was performed on 200 ng of gDNA using the EZ DNA Methylation-Lightning Kit (Zymo Research, Irvine, CA, USA), following the manufacturer’s instructions. The resulting bisulfite-converted DNA (bcDNA) was eluted in 15 μl of buffer, yielding 13.33 ng/μl, assuming 100% recovery [57]. Sensitivity was assessed by serial dilutions down to 12.5 ng of gDNA input, resulting in 0.83 ng/μl of bcDNA. To validate methylation accuracy, standards were created by mixing fully methylated (100%) and unmethylated (0%) DNA in 10% increments using the EpiTect PCR Control DNA Set (Qiagen).

Primer design for methylation analysis

Bisulfite sequencing primers were designed using Pyromark Assay Design Software version 2.0 (Qiagen) and Bisearch (http://bisearch.enzim.hu), with optimal parameters including a melting temperature range of 56–62°C, 100 bp amplicon size, 25 bp primer length, and exclusion of degeneracy and secondary structures (Supplementary Table 6). SBE primers with extended T-tails for multiplex analysis were designed (Supplementary Table 7), and the final primer sets used in the model are detailed in Supplementary Table 8.

Bisulfite sequencing using single base extension

For the SBE assay, 1 μL of bcDNA (13.33 ng/μL) was amplified in a 20 μL multiplex PCR reaction containing 2 μL of 10× Gold ST*R Buffer (Promega, Madison, WI, USA), 3.5 U of AmpliTaq Gold DNA Polymerase (Applied Biosystems, Foster City, CA, USA), and sequence-specific primers. The PCR protocol began with an initial denaturation at 95°C for 11 minutes, followed by 34 cycles of 94°C for 20 s, 56°C for 60 s, and 72°C for 30 s, concluding with a final extension at 72°C for 7 min. The amplified PCR products were first purified using EXOSAP-IT^™ (Applied Biosystems) at 37°C for 45 min, then heat-inactivated at 65°C for 15 min. DNA methylation detection by SBE was performed in a 10 ul reaction mixture containing 1 μL of 10× SNaPshot Multiplex Ready Reaction Mix, 2 μL of 5× BigDye^™ Terminator v.1.1 and 3.1 Sequencing Buffer (both from Applied Biosystems), SBE primers, and 1 μL of purified PCR product. SBE thermal cycling conditions consisted of 96°C for 10 s, 50°C for 5 s, and 60°C for 30 s, repeated for a total of 28 cycles. The SBE products were then purified using recombinant SAP (Applied Biosystems) under the identical conditions as the previous purification. All PCR and SBE reactions were conducted on a Veriti 96-Well Thermal Cycler (Applied Biosystems, Waltham, MA, USA), and sequencing was performed using a 3500 Genetic Analyzer (Thermo Fisher Scientific) equipped with a 36 cm capillary and POP-4 polymer.

Bisulfite sequencing using next generation sequencing

For targeted NGS, libraries were prepared using the KAPA HyperPrep Kit, KAPA Universal Adapters, and KAPA Unique-Dual Indexed (UDI) Primer Mixes, in accordance with the manufacturer’s guidelines provided in the KAPA HyperPrep Kit technical datasheet. Initially, 13.33 ng of bcDNA underwent targeted multiplex PCR with specific primer sets (Supplementary Table 3). PCR products exceeding 300 ng were then subjected to end-repair and A-tailing procedures. The reaction mixture was incubated with a universal adapter at 20°C for 1 hour. Following adapter ligation, each library was purified using AMPure XP beads (Beckman Coulter, Brea, CA, USA). Libraries were then amplified using UDI primer mixes through 13 cycles of PCR and subsequently purified again with AMPure XP beads. The concentration and quality of the libraries were assessed using the 2100 Bioanalyzer and 4150 TapeStation (Agilent, Santa Clara, CA, USA), and all samples were pooled into a single batch for further processing. Double-size selection was performed using SPRIselect beads (Beckman Coulter) to ensure the desired fragment size. The final library was quantified using KAPA Library Quantification Kits (Roche, Basel, Switzerland), following the provided instructions. The pooled library was adjusted to 4 nM, denatured with 0.2 N NaOH, and adjusted to 20 pM with Illumina prechilled hybridization buffer. The denatured library was subsequently diluted to 7–9 pM with a spike-in of 2.5% PhiX control v.3. Sequencing was performed on an Illumina MiSeq system using the MiSeq Reagent Kit v.3, with a 2 × 300 cycle configuration (Verogen, San Diego, CA, USA).

Methylation data processing

SBE methylation data were analyzed using GeneMapper Software 5 (Thermo Fisher Scientific). Methylation values were calculated from the peak height, representing the relative fluorescence units (RFU) of C or G nucleotides corresponding to methylated DNA, normalized by the total DNA, determined as the sum of C+T or G+A nucleotides. Raw sequencing data from NGS (FASTQ files) were processed on CLC Genomics Workbench 21.0.5 (Qiagen). The analysis workflow included merging paired-end reads, trimming low-quality bases based on Phred scores reflecting a 1% error probability, and excluding reads containing two or more ambiguous bases or measuring less than 10 bp in length. Processed reads were aligned to the reference genome (Human GRCh38), and methylation calling was performed, omitting loci with coverage below 20. The resulting methylation values were utilized for subsequent machine learning modeling.

Machine learning

Machine learning analyses to develop age prediction models were performed using Python 3.8.8. Linear regression (LR), stepwise linear regression (SLR), Lasso Regression (Las), ridge regression (Rid), and elastic-net regression (Ela) were implemented. Tree-based decision models, including XGBoost (XGB), gradient boosting (GB), and random forest (RF), were implemented using the Scikit-learn and XGBoost libraries in Python. Hyperparameters were optimized with GridSearchCV, and model performance was evaluated with metrics including the Pearson correlation coefficient (r), R-squared (R²), mean absolute error (MAE) and root mean square error (RMSE).

Statistical analysis

Statistical comparisons for mRNA expression were performed using the Chi-square test and the unpaired t-test. Data are presented as mean ± standard error of the mean (SEM) and were analyzed using GraphPad Prism 5.01 (GraphPad, La Jolla, CA, USA). Levels of statistical significance are denoted by asterisks with the following thresholds: ^*p ≤ 0.05; ^**p ≤ 0.01; ^***p ≤ 0.001.

Availability of data

All data supporting the findings of this study are available within the paper and its Supplementary files. Raw data underlying the main and supplementary figures and tables have been deposited in the Gene Expression Omnibus (GEO) under accession numbers GSE244996 and GSE294234. A source data file containing statistical analyses and p-values is provided as Supplementary Table 1.