Introduction

The field of genetics focuses on heritable alterations in gene function or activity that result from direct modifications to the DNA sequence. These changes include translocations, insertions, deletions, and point mutations. Conversely, the field of epigenetics investigates heritable alterations to the function or activity of genes that transpire in the absence of modifications to the fundamental DNA sequence1,2. Despite having the same genetic content in any given organism, different cell types express different genes at different times. More broadly, in multicellular organisms, diverse gene expression is mediated in a range of cells and tissues by several factors, including epigenetic mechanisms such as DNA methylation, histone modifications, chromatin remodeling, and non-coding RNA-mediated regulation3.

One of the most studied and mechanistically understood epigenetic modification is the DNA methylation, which is crucial to the development of mammals4, and is conserved in the majority of plants, animals, and fungi5. DNA methylation is regulated by DNA methyltransferases, methyl-CpG binding proteins, and other chromatin-remodeling agents6. It is described as the addition of a methyl (CH3) group into DNA, and it mostly happens on CpG sites or cytosine that precede guanine nucleotides in mammalian cells7. The impact of DNA methylation can vary based on their ___location in the genome. Methylation at promoter regions typically suppresses gene expression, while enhancing it when found within the gene sequence8. DNA methylation changes chromatin structure and other epigenetic markers, like histone modifications and transcription factors, thus controlling gene expression8,9. An increasing prevalence of human diseases is linked to improper establishment and/or maintenance of this epigenetic information, emphasizing the significance of DNA methylation10,11. Moreover, there is growing interest in creating pharmacological interventions to correct epigenetic abnormalities12.

The importance of comprehensively studying the human DNA methylation in relation to human health and disease has been acknowledged with the advancement of technologies that enable the measurement of DNA methylation at both the individual and population levels across significant regions of the human genome. High-throughput technologies such as the Illumina Infinium HumanMethylation450 BeadChip13,14, the MethylationEPIC BeadChip7, allow for the assessment of methylation at hundreds of thousands of CpG sites across the genome. Additionally, next-generation sequencing (NGS) techniques, including whole-genome bisulfite sequencing (WGBS)15 and reduced representation bisulfite sequencing (RRBS)16, provide detailed and comprehensive methylation profiles. In this study we used the newly developed “Infinium HumanMethylationEPIC version 2.0”17, an array including over 920,000 CpG probes.

In the present study, we investigate whether different isoforms of a growth hormone receptor (GHR) gene can exhibit distinct DNA methylation patterns and contribute to epigenetic variation. The GHR gene can appear in two isoforms: the full-length wild-type (WT) GHR (flGHR), and a shorter gene with a deletion of exon 3 (d3GHR); this common polymorphic deletion variant was described in 2000, and is associated with increased receptor activity as a result of improved signal transduction18. The distribution of GHR isoforms includes, on average, 50% of the population with flGHR (fl/fl), 30–40% that are heterozygous for the exon 3 deletion (fl/d3), and 10–20% that are homozygous for the deletion (d3/d3)19,20,21,22. According to our previous study23, exon 3 deletion is associated with increased height, lower serum IGF-1 levels, and enhanced GH sensitivity. Moreover, the d3/d3 genotype extends life expectancy in males by about 10 years, based on multivariate regression analysis23. Building upon these fascinating discoveries, we carried out an in-depth analysis of DNA methylation in order to gain a better understanding of this isoform’s impact on aging, with the anticipation of uncovering differences in aging trends and biological age estimates due to genetic differences. Furthermore, according to published studies, this deletion was associated with more gradual age-related weight and height gains24, greater spontaneous growth25,26, and a better response to human growth hormone treatments in a variety of growth disorders, compared to individuals with the full-length genotype19,27,28.

This project involves individuals from different Druze clans residing in various distinct villages in the northern region of Israel. The Druze religion has its roots in the eleventh century, and they are known for being a closed and endogamous population with a high level of consanguinity, as the common practice is first cousins marriages29. Isolated genetic populations offer a valuable resource for research focused on the genetic predisposition behind diseases and their characteristics.

In the current study, we studied the DNA methylation profiles of individuals with different d3GHR isoforms. Given that DNA methylation alterations are essential for regulating gene expression and body functions30, we hypothesized that individuals with the homozygous deletion will show differentially methylated positions (DMPs) and regions (DMRs) compared with the heterozygous and WT groups, and will have a beneficial DNA methylation pattern which promotes positive gene expression regulation, cellular differentiation, and overall genomic stability, and exhibit younger biological age. Our findings indicate that varying levels of methylation are found among individuals based on their genotype, and that this is associated with exceptional biological and chronological age. To conclude, this study contributes to the growing field of epigenetics and provide insight into how genetics and epigenetic mechanisms impact biological diversity and health outcomes in individuals carrying at least one d3 allele of GHR.

Results

Subject characteristics and genotyping

We analyzed the genotypes of 73 individuals from six families living in four Druze villages/townships, notably, these families were unrelated to one another, ensuring diverse genetic backgrounds for our investigation. The findings indicated that 43 individuals (46.57%) possessed the WT GHR genotype (fl/fl), 33 individuals (45.2%) were carriers (heterozygous) of the exon 3 deletion (d3/fl), and 6 individuals (8.21%) had two copies (homozygous) of the exon 3 deletion (d3/d3). The maximum age observed in the WT group was 104 years, with the youngest participant being 23 years old; in the heterozygous group, ages ranged from 21 to 100 years; and in the homozygous group, participants ranged in age from 43 to 109 years. These age ranges provided a diverse representation of age groups within each genotype, enabling a comprehensive analysis of age-related differences in DNA methylation patterns. The d3GHR genotype is more common in older individuals than in younger individuals, and its prevalence increases with age (P < 0.05), with males showing an increase from 9 (ages 41–65) to 10% (ages 66–85), and a substantial rise to 20% in the 90 + age group. Similarly, females exhibit an increase in d3GHR frequency from 7% (ages 41–65) to 16.6% (ages 90 +), indicating pronounced age-related variations in genotype distribution among both genders (Fig. 1).

Fig. 1
figure 1

Prevalence of d3GHR homozygotes among (a) males (n = 40), and (b) females (n = 33). *P < 0.05.

Full size image

A group of 14 individuals was chosen for additional examination of DNA methylation patterns. This group included individuals with WT (n = 4), heterozygous (n = 6), and homozygous (n = 4) genotypes (Supplementary Table 1), belonging to six distinct families (Fig. 2). The 14 individuals chosen for DNA methylation analysis were selected based on factors, such as sample availability, family relationships, genetic diversity, and representation of target genotypes to thoroughly study epigenetic variations associated with the investigated genotypes. The familial lineage of these individuals encompassed information on ages and genotypes for all collected samples. This enabled us to analyze and contrast DNA methylation patterns among various genotypic backgrounds.

Fig. 2
figure 2

Representation of the familial relationships, ages, and genotypes of the individuals included in our study among six separate clans (a–f). Black arrows indicate the individuals who were included in the DNA methylation analysis, a diagonal line indicates individuals who have died (Abbreviations: YO—years old).

Full size image

DNA methylation

DNA methylation data from all 14 individuals (marked in Fig. 2 and in supplementary Table 1) underwent quantile normalization. Quantile normalization normalizes the methylation beta values (β_meth) across samples, therefore all samples get similar β_meth values (Fig. 3A). Mean β_meth values were calculated for each sample based on beta values of CpG sites that passed QC filtering, serving as an indicator of overall DNA methylation levels (Fig. 3).

Fig. 3
figure 3

Average methylation level comparison across groups. Methylation levels were analyzed using the Infinium MethylationEPIC Array and converted to beta values using preprocessQuantile to eliminate background read signal with Minfi R package. (a) Average β_meth values per group, WT = 0.5778, Heterozygous = 0.5788, and Homozygous = 0.5781, no significant difference between groups were identified. (b) Linear regression analysis of average β_meth values across age (p < 0.05).

Full size image

Differential methylation analyses- DMP and DMR

We analyzed the three separate groups by looking at their DNA methylation patterns in order to identify potential links to health conditions. Examination of methylation data across the genome showed significant variations in DNA methylation patterns between the groups, especially in CpG sites linked to biological functions. To examine the epigenetic factors linked to exon 3 deletion, we compared methylation levels in individuals with homozygous flGHR, heterozygotes, and homozygotes for the exon 3 deletion. We identified five unique methylated position DMPs which were notably significant (Table 1) (P-values have been adjusted using the Benjamini–Hochberg procedure, and are labeled as ‘adj.p-value (BH)). In this analysis, we only included male subjects and combined heterozygous and WT individuals (since there is only one male with the WT genotype) to compare with d3 homozygous individuals. The unique patterns seen indicate potential epigenetic signatures specific to each group. In the d3GHR group, three DMPs exhibited hypermethylation: cg09306641 on chr18, cg09508496 on chr19, and cg23837191 on chr18. Furthermore, two DMPs showed hypomethylation compared to the WT/Heterozygous group: cg26168992 on chr5 and ch.14.569218R on chr14 (Padj < 0.05) (Fig. 4).

Table 1 Top 5 significant DMPs that were identified.
Full size table
Fig. 4
figure 4

5 Differentially Methylated Positions (DMPs) that were identified. Displaying the DMP value for each group with hyper- or hypo-methylation. The x-axis indicates the GHR isoform genotype and the y-axis represents the methylation beta value of each CpG site. The WT/heterozygous and homozygous groups are depicted by red and blue boxes, respectively. *p < 0.05, **p < 0.005.

Full size image

Next, we proceeded with the identification of differentially methylated regions (DMRs) and found two significant DMRs (Table 2), each composed of two and three CpG sites, respectively (FDR < 0.05) (Table 3 and Fig. 5). Two CpGs are located in proximity to each other on chromosome 18 (cg23837191 and cg09306641), while another three are located in close proximity on chromosome 19 (cg02452418, cg02452966 and cg09508496). The DMR distributions and their significance to the chromosomes are shown in Table 2. The homozygous d3/d3 group had higher methylation levels than the WT and heterozygous groups in all CpG sites (Fig. 5). Our findings revealed two DMRs located on chromosome 18 and demonstrated a hypermethylation pattern. These DMRs are adjacent to the Nuclear Factor of Activated T Cells 1 (NFATC1) gene (chr18: 79,395,930-79,529,323), which is involved in the activation of cytokine gene expression in T-cells, particularly in driving IL-2 or IL-4 gene transcription. It also regulates gene activation in developing cardiomyocytes31. Additionally, another DMR locates next to the RYR1 gene (chr19: 38,433,69-38,587,564), which plays a crucial role in muscle function and regulation of calcium signaling32. This DMR is located on chromosome 19 at genomic coordinates 38,565,444-38,565,576, and exhibits hypermethylation in our analysis compared to WT/Heterozygous samples.

Table 2 Top 2 significant DMRs identified between the GHR isoforms.
Full size table
Table 3 DMPs identified in the two significant DMRs.
Full size table
Fig. 5
figure 5

DNA methylation levels of 2 differentially methylated regions (DMRs). Five genomic positions were identified. The x-axis indicating the GHR isoform genotypes and the y-axis represents the adjusted methylation beta value of each CpG site. The WT/heterozygous and homozygous groups are depicted by red and blue boxes, respectively. *p < 0.05, **p < 0.005.

Full size image

These findings highlight distinct epigenetic signatures associated with each group, providing insights into the underlying molecular mechanisms that may contribute to their unique physiological and pathological characteristics.

To visualize the locations of the DMRs identified in our study, we constructed a genetic map of the relevant chromosomes. This map highlights the precise positions of each DMR, providing a clear overview of the methylation changes across the chromosomal landscape (Fig. 6). In specificity, the genetic map (Fig. 6b) shows the exact ___location of DMR2 in relation to the RYR1 gene.

Fig. 6
figure 6

The plot illustrates the genetic map of the specified chromosomes, highlighting the positions of the identified Differentially Methylated Regions (DMRs) identified by the DMRcate analysis, the positions of nearby genes (including DMR2 in relation to the RYR1 gene), and the methylation levels of individual samples displayed as a heatmap. DMRs are indicated by violet bars and marked with black arrows, while the RYR1 gene is represented by orange bars. Genomic coordinates are based on data obtained from the UCSC Genome Browser.

Full size image

DNAm age clocks

Epigenetic clocks allow comparison of an individual’s chronological age to their biological age, as calculated based on DNAm clock CpGs. Estimators of age based on DNA methylation, include the pan-tissue epigenetic clock by Horvath 2013, which is based on the values of 353 distinct CpG sites to create a multivariate age predictor33, and an estimator that was created based on 71 CpGs in leukocytes by Hannum 201334. Following Horvath’s first DNA methylation clock ‘HorvathAge’, additional second-generation methylation clocks have been created, e.g. GrimAge and PhenoAge, to focus on CpG sites that may play a more significant role in aging. All age clocks aim to estimate an individual’s biological age based on molecular markers and offer insights into an individual’s biological aging processes and potential health outcomes. This study uncovered an interesting finding about the connection between biological age and chronological age among the study participants. Our findings consistently showed that the predicted biological age was lower than the real chronological age for the majority of participants studied, uncovering intriguing patterns within the genotypic groups that were analyzed.

According to our results, the mFitAge biological age clock revealed significant differences between chronological age and biological age among all individuals assessed (P < 0.05). Notably, as individuals progressed in chronological age, the disparity from biological age increased, indicating a widening gap between actual age and the biological aging processes. This trend was consistent across all biological age clocks utilized in our analysis, emphasizing the robustness of the observed phenomenon.

The “DNAmFitAge” calculator indicated varied biological ages in all genotypic groups (Supplementary Table 2). More precisely, the WT and heterozygous group had a mean biological age of 93.968 years, while the homozygous group had a mean biological age of 89.104 years. The delta age, which represents the difference between chronological and biological age, revealed a notable trend favoring the homozygous d3GHR group. The results revealed a positive delta age of + 4.229 years for the d3GHR group, indicating a possible deceleration or more favorable aging phenotype in which biological age was lower than chronological age. In contrast, the heterozygous and WT groups had negative delta ages of -0.593 years. The observed delta age differences were statistically significant (P < 0.05), indicating distinct aging trajectories associated with different genotypic backgrounds (Fig. 7). The observed increase in delta age among the d3GHR group highlights the potential influence of genetic variations, such as the exon 3 deletion in the growth hormone receptor (GHR) gene, on aging mechanisms and estimations of biological age. This finding underscores the importance of considering genetic differences in assessing biological age and emphasize the importance of considering both chronological age and biological age metrics in aging-related research and health outcomes.

Fig. 7
figure 7

Comparison of the delta age (chronological age and biological age) among various GHR genotypes of individuals aged 65 and above. The x-axis represents the different genotypes (WT and heterozygous vs. homozygous d3), the y-axis represents the delta age values, calculated as the discrepancy between a person’s chronological age and biological age (p < 0.05).

Full size image

To conclude, the discovery of a consistently younger biological age brings up fascinating inquiries about the various factors that could be causing this occurrence, such as genetics, phenotypes, environment, and lifestyle choices.

Discussion

Epigenetics is becoming increasingly common among researchers in various fields, including cancer, molecular medicine, behavior, development, physiology, morphology, and cell biology35,36,37,38. This research field plays a pivotal role in regulating cellular function, and encompasses a wide range of research areas focused on understanding how gene expression and phenotypic characteristics are regulated and influenced by factors beyond the DNA sequence itself39. DNA methylation, the most abundant epigenetic regulator, is crucial for controlling genes and differentiating cells, and has become a field of interest in complex disease processes due to its ability to either direct or mirror transcriptional activity40,41. Variability in DNA methylation levels in healthy individuals is also suggested to impact susceptibility to common diseases and response to therapies42,43. Furthermore, variations in DNA methylation have also been observed among different ethnic groups44,45.

In a previous study, we found that men who are homozygous for the d3GHR isoform live 10 years longer on average23. Additionally, according to published studies, the d3GHR has a positive impact on growth disorders, including increased spontaneous growth19,25,46,47. Building on these fascinating results, we carried out a thorough DNA methylation analysis to learn more about this isoform’s involvement in the aging process. Our research sought to pinpoint the d3GHR isoform’s differentially methylated positions (DMPs) and regions (DMRs), in order to shed light on the epigenetic processes that underlie its lifespan-extending benefits. Furthermore, we compared the biological age of homozygotes for the d3GHR isoform to those with the full-length gene, hoping to uncover potential differences in aging trajectories and biological age estimates caused by these genotypic variations in a family-based structure. This study is a step toward understanding the complex interplay between genetic variations, epigenetic modifications, and aging-related phenotypes, which have implications for longevity-related mechanisms and age-related diseases.

Here we report on methylation profile analyses comparing between individuals homozygous for the full length GHR (fl/fl) and heterozygotes, and individuals homozygous for the exon 3 deletion isoform (d3/d3). The participants in our study are from the Druze community. The Druze community in Israel is a distinctive religious and ethnic minority found mainly in the northern part of the country. Due to their high rate of consanguinity, it is more common to observe recessive traits, which allows the family to maintain specific isoforms or traits across generations. Consequently, in this study, if the oldest individual (proband) in such a lineage possesses the d3GHR isoform, this isoform is likely to be observed in subsequent generations, including children and grandchildren, in both homozygous and heterozygous forms. This inheritance pattern highlights the impact of familial genetic relationships on the preservation of particular genetic traits and ensuring the preservation and transmission of such recessive traits within the family over successive generations.

Among 73 Druze clan members in this research, 8.2% had the d3GHR isoform, with almost 11% being males. A noticeable age-related rise (p = 0.04) in this isoform was observed among males. In our recently published research on Ashkenazi Jews, a noticeable variation in allele frequency for the d3GHR polymorphism was observed in Ashkenazi males. While the male control group had only 4% homozygote deletions, male centenarians and their offspring carried 12 and 11%, respectively23.

Our analysis revealed a significant decrease in DNA methylation levels with advancing age across the study cohort (Fig. 3B). This age-related decrease in methylation was observed at specific CpG sites, as well as in broader genomic regions, indicating a widespread epigenetic remodeling process associated with aging. Our results align with and are backed up by many studies showing a decline in DNA methylation in different populations and tissues as individuals grow older48,49,50. Moreover, after reaching adulthood, numerous studies have observed an average reduction in blood DNA methylation as individuals grow older34,51,52. This reduction can result in losing control over transcription and can lead to or contribute to negative effects of aging53.

In the present study, despite the constraints of a small sample size, we found five significant DMPs with unique methylation patterns. These DMPs exhibited reliable variations in DNA methylation levels, underscoring their possible biological significance. These discoveries set the foundation for future research on the functional significance of these DMPs and their possible involvement in disease development or biological functions. Furthermore, our analysis revealed the presence of two significant DMRs characterized by distinct methylation patterns. Notably, one DMR consisted of two CpG sites located on chromosome 18, while the other encompassed three CpG sites located on chromosome 19, each showing significant differences in DNA methylation levels in the d3 homozygous samples compared to control samples (WT and heterozygotes).

We identified two DMRs, which are constituted by five DMPs (Table 3), those DMPs are located adjacent to NFATC1 and RYR1 genes. Notably, these positions showed hypermethylation patterns among the d3 homozygous group. This finding suggests a regulatory association between DNA methylation status and the genetic loci that surround NFATC1 and RYR1. The observed hypermethylation near these genes in d3 homozygous individuals highlights potential epigenetic mechanisms that may influence gene expression and cellular processes related to NFATC1 and RYR1, necessitating additional research into their functional implications and role in biological pathways.

Regarding the DMPs that were identified in Table 1, although their functions are unknown, the notable differences in methylation at these sites indicate their possible significance in epigenetic control and cellular functions. Continued research efforts to understand the importance of these non-coding or unannotated regions could reveal new knowledge about the intricacies of the epigenome and its role in different conditions.

As recognition of the important impact of epigenetic changes on the aging process increases, DNA methylation patterns are now used to determine biological age, known as the epigenetic clock. Levels of DNAm have been used to create precise composite biomarkers for measuring biological age33,34,54,55. Studies of epigenetic age acceleration in age-related disorders have been conducted using various epigenetic clocks. Multiple studies indicated that DNA methylation age acceleration is linked to the occurrence, future development, and death rates of various cancer types56,57,58. Individuals may have a biological age that is younger than their actual chronological age, since age-associated methylation changes can affect the expression of genes33, and distinct DNA methylation profiles may contribute to extended lifespan, these profiles often include the maintenance of youthful methylation patterns in genes related to immune response, inflammation, and metabolic regulation59. In our study, the DNAmFitAge clock has yielded compelling results regarding the delta age, representing the difference between chronological age and biological age, which consistently favored the homozygous group. This finding aligns perfectly with our initial hypothesis, suggesting that genetic variations, specifically the homozygous d3GHR state, may contribute to a more favorable aging phenotype characterized by a lower biological age compared to chronological age. This underscores the potential impact of genetic factors, such as the presence of specific alleles or genetic variants, on aging processes and biological age estimates. The observed positive trend in delta age provides further support for the hypothesis that the exon 3 deletion of the GHR gene, represented by the homozygous group in our study, may exert beneficial effects on aging trajectories.

In conclusion, the inclusion of a family structure sample size enhances the reliability and generalizability of our findings, allowing for more comprehensive insights into DNA methylation patterns and their associations with various biological processes and phenotypes. Our results highlight the significance of studying different genotypes in such structure, as they may lead to various genetic effects, particularly in influencing DNA methylation patterns. Furthermore, by considering genotype-specific analyses, we can gain deeper insights into the interplay between genetic variations and epigenetic mechanisms, providing a more nuanced understanding of how genetic factors contribute to complex biological processes and disease susceptibilities. These considerations are crucial for advancing precision medicine approaches and developing targeted interventions aimed at modulating epigenetic regulation for improved health outcomes.

Materials and methods

Sample collection and preparation

We recruited long-lived Druze probands (age > 90 years), as well as available spouses, children and grandchildren. In accordance with the principles outlined in the Helsinki committee approval, all procedures involving human participants were approved by “Rambam Health Care Campus, Haifa, Israel, 0038-14-RMB”. Written informed consent was obtained from all participants prior to their inclusion in the study. Family history information, including ages, morbidity and mortality, was collected in detail from each family, and all participants were evaluated by medical professionals. We collected 9 mL of EDTA whole blood for genomic DNA (gDNA) extraction, using the FlexiGene DNA Kit (Qiagen)DNA concentration was quantified using the DS-11 spectrophotometer (DeNovix).

Multiplex PCR and samples genotyping

Individual genotypes were assessed for the d3GHR polymorphism. A pair of primers were used to amplify the extracted DNA via PCR. The G1 and G2 primers, which encompass exon 3, result in a PCR fragment of 532 base pairs for the deletion (d3) allele, while the G1 and G3 primers yield a PCR fragment of 935 base pairs for the full-length WT (fl) allele (G3 is an antisense primer located inside exon 3) (GenBank accession number AF155912) (Supplementary Table 3). The PCR reaction mixture contained 0.5 µl of forward primer and 0.5 µl of reverse primer, 10 µl Taq X2 Master Mix, 1 µl template DNA, and 8 µl nuclease-free water. Utilizing a Bio-Rad device, the following parameters were used for the PCR process: 94 °C for 5 min, 35 cycles of 94 °C for 30 s, 60 °C for 30 s, 72 °C for 90 s, and a final extension at 72 °C for 7 min. A 1.5% ethidium bromide stained agarose gel stained was used to separate the PCR products by size (Supplementary Fig. 1). To successfully confirm the deletion, the 953 and 532 bands were cut and purified using a GeneJET Gel Extraction Kit (Thermo Fisher Scientific, Cat. No. K0692), followed by Sanger sequencing by Macrogene sequencing services.

EPIC-methylation assay and bioinformatics analysis

Extracted gDNA was sent to Diagenode for epigenome DNA methylation analysis using the Infinium MethylationEPIC Array Service v2.0. This array is an enhanced version of EPIC v1.0, the new cutting-edge DNA methylation analysis technique by Illumina that utilizes bisulfite conversion. It permits the quantitative measurement of the overall methylation status of more than 930,000 methylation sites across the human genome with a resolution of single nucleotide. It provides an extensive range of information on CpG islands, enhancer regions, open chromatin sites, and other key areas of the methylome (Diagenode Cat# G0209006). Briefly, data analysis was performed within the R v4.2.2 statistical environment, raw intensity (.idat) files were imported into R using the minfi package60, followed by the calculation of M and methylation Beta values. Methylation beta (β_meth) values refer to the proportion of methylation at a specific CpG site, ranging from 0 (no methylation) to 1 (complete methylation). A detection p-value was computed to assess signal quality at specific sites. Stratified quantile normalization of beta and M values was performed using the preprocessQuantile function in the minfi package. Density plots were generated to assess the distribution of beta values before and after quantile normalization (Supplementary Fig. 2). Only male participants were considered for the differential methylation examination (n = 9), and the comparison between “homozygous versus heterozygous/wild-type” was employed to increase statistical strength in all models examining differential methylation (DMP and DMR). A flow chart diagram depicting the analysis of the “Infinium MethylationEPIC Array” is presented in Supplementary Fig. 3.

Differentially methylated probes (DMPs)

DMPs associated with a genotype were identified using M values and a linear model from the limma R package (v.3.54.2)61. The effect of the ‘genotype’ coefficient was extracted to identify significant DMPs, with all DMPs having an adjusted p-value under 0.05 considered significant, lmFit function from limma R package was used. Adjustment was performed with Benjamini–Hochberg procedure. The Benjamini–Hochberg adjustment is a statistical method that controls for false discovery rate (FDR).

Differentially methylated regions (DMRs)

For the identification of DMRs spanning multiple CpG sites, the DMRcate R package (v.2.12.0) was employed62. All DMRs that met the criterion of a false discovery rate (FDR) below 0.05 were deemed statistically significant.

Probes were annotated using the IlluminaHumanMethylationEPICv2anno.20a1.hg38 package [https://github.com/jokergoo/IlluminaHumanMethylationEPICv2anno.20a1.hg38], distinguishing between Type I and Type II probes and annotating them to CpG Islands, open sea regions, shores and shelves.

Creation of genetic map

To illustrate the locations of the Differentially Methylated Regions (DMRs) in their genomic context, we used the “DMR.plot” function. This function allows for the creation of detailed genetic maps that display the positions of DMRs relative to the chromosome and to nearby genes. The resulting plots provide a clear visual representation of how these DMRs are distributed along the genome and their proximity to key genetic features.

Epigenetic clock and DNA methylation (DNAm) age calculation

We utilized epigenome-wide methylation information to calculate the epigenetic clocks using the DNA Methylation Age Calculator (https://dnamage.genetics.ucla.edu). The data used for the calculation were adjusted using the DNA methylation Age Calculator.