Epigenome-wide association study identifies DNA methylation loci associated with handgrip strength in Chinese monozygotic twins

Background: The decline in muscle strength and function with aging is well recognized, but remains poorly characterized at the molecular level. Here, we report the epigenetic relationship between genome-wide DNA methylation and handgrip strength (HGS) among Chinese monozygotic (MZ) twins. Methods: DNA methylation (DNAm) profiling was conducted in whole blood samples through Reduced Representation Bisulfite Sequencing method. Generalized estimating equation was applied to regress the DNAm of each CpG with HGS. The Genomic Regions Enrichment of Annotations Tool was used to perform enrichment analysis. Differentially methylated regions (DMRs) were detected using comb-p. Causal inference was performed using Inference about Causation through Examination of Familial Confounding method. Finally, we validated candidate CpGs in community residents. Results: We identified 25 CpGs reaching genome-wide significance level. These CpGs located in 9 genes, especially FBLN1, RXRA, and ABHD14B. Many enriched terms highlighted calcium channels, neuromuscular junctions, and skeletal muscle organ development. We identified 21 DMRs of HGS, with several DMRs within FBLN1, SLC30A8, CST3, and SOCS3. Causal inference indicated that the DNAm of 16 top CpGs within FBLN1, RXRA, ABHD14B, MFSD6, and TYW1B might influence HGS, while HGS influenced DNAm at two CpGs within FBLN1 and RXRA. In validation analysis, methylation levels of six CpGs mapped to FLBN1 and one CpG mapped to ABHD14B were negatively associated with HGS weakness in community population. Conclusion: Our study identified multiple DNAm variants potentially related to HGS, especially CpGs within FBLN1 and ABHD14B. These findings provide new clues to the epigenetic modification underlying muscle strength decline.


Introduction
Muscle strength is a crucial factor in healthy aging and an important predictor of several adverse health outcomes, which decline from midlife and throughout life (McLeod et al., 2016;Balogun et al., 2019;Alcazar et al., 2020).HGS, a simple measure of upper body muscle strength in clinical practice (Lauretani et al., 2003), is influenced by both environmental and genetic factors (Tian et al., 2017).The correlation between HGS and various factors, encompassing nutritional parameters and levels of physical activity, has been extensively explored through epidemiological studies (Rondanelli et al., 2016;Mendes et al., 2017;Weng et al., 2022).Moreover, previous genome-wide association studies (GWAS) have elucidated the influence of genetics on HGS, revealing heritability estimates within the range of 40%-65% (Tiainen et al., 2004;Matteini et al., 2010;Matteini et al., 2016).Recent emerging GWASs have facilitated the identification of human genetic variants associated with HGS (Matteini et al., 2016;Luo et al., 2022).Nevertheless, the single nucleotide-level polymorphisms (SNPs) reported in these studies only explained a small portion of heritability of HGS, suggesting that additional gene-regulatory mechanisms are involved in muscle maintenance or decline (He et al., 2020).
Emerging evidence indicated that epigenetic processes play a crucial role in the development of complex human traits or diseases (Godfrey et al., 2007).Epigenetic modifications serve as potential mediators for the impact of external factors on muscle maintenance or decline, providing a mechanistic framework that integrates environmental and genetic influences.DNA methylation, a wellstudied epigenetic modification, provides a novel perspective on the variability of muscle strength.It is widely recognized that DNA methylation is a fundamental mechanism underlying myogenesis (Carrió et al., 2015;Laker and Ryall, 2016).Previous research has suggested that DNAm regulates skeletal muscle differentiation and proliferation by modulating the expression of genes from the myocyte enhancer and myogenic pathway (Bharathy et al., 2013).Several studies have reported a hypermethylated state in skeletal muscle tissue in older subjects compared to younger subjects (Zykovich et al., 2014;Turner et al., 2020).There are multiple epigenome-wide association studies (EWAS) of sarcopenia or muscle mass and strength among people of European ancestry (Livshits et al., 2016;Voisin et al., 2021;Antoun et al., 2022).These studies have identified many CpG loci and differentially methylated regions (DMRs) that may be associated with muscle strength.Nevertheless, three EWASs have tested associations between whole-blood DNAm level and HGS, but no significant association between methylation levels at single CpG and HGS was observed (Bell et al., 2012;Marioni et al., 2015;Soerensen et al., 2019).
The trait-discordant monozygotic twins design is a valuable and powerful design for EWAS that can link epigenetic modification to traits, as the genetic component is perfectly matched within monozygotic twin pairs (Tan et al., 2015).This design also provides an opportunity to identify different DNA methylation patterns triggered by environmental factors.Due to difference in environmental exposures, genetic backgrounds, and lifestyles among different ethnicities, DNAm patterns are likely to exhibit distinctions.Hence, it is necessary to conduct comparable studies using samples from different ethnic groups.Currently, relevant investigation has not been reported in the Chinese population.

Participants and study procedures
The study population was a subsample of the Qingdao Twin Registry System.For detailed information on recruitment procedures, please refer to our previous work (Xu et al., 2017).Participants were excluded if they suffered from stroke, cardiovascular disease, and/or tumor, or were unable to complete the examination.Twin pairs with intra-pair HGS difference ≥0.1 kg were chosen for the trait-discordant monozygotic twin design based on our previous experience with similar studies (Wang W. et al., 2021;Wang T. et al., 2021;Wang et al., 2023).A total of 66 HGSdiscordant monozygotic twin pairs were included.The Regional Ethics Committee of the Institutional Review Committee of Qingdao CDC approved this study.This study followed the Helsinki Declaration and all participants provided informed consent.
After a fasting period of 10-12 h, the participants completed a questionnaire and physical examination, and a venous blood sample of 10 mL was collected.Handgrip strength was measured using a hand-held dynamometer (WCS-100, Nantong, China).Participants were asked to forcefully squeeze the dynamometer three times for each hand, with the maximum value being used for subsequent analysis.

DNA methylation data
Total DNA was extracted from venous blood samples for reduced representation bisulfite sequencing experiments.In summary, the initial step involved digestion of genomic DNA to produce shorter fragments.Subsequently, the CpG-rich DNA fragments underwent bisulfite conversion.Finally, the resulting cDNA library was sequenced.The resulting raw methylation data encompassed 551,447 CpGs throughout the genome of each individual.The raw reads were mapped to the human Genome Reference Consortium Human Build 37 using Bismark (Krueger and Andrews, 2011), and methylation levels were determined by smoothing the data with the BiSeq package in R software.We maintained control over the coverage, ensuring it was within the 90% quantile.CpGs with an average methylation β-value below 0.01 or those with more than 10 missing observations were excluded.The methylation β-values were transformed to M-values by log 2 transformation (Wang T. et al., 2021).

Cell-type composition estimate
DNAm patterns vary across different cell types.To address the potential confounding effect of cell-type composition on DNAm analysis in whole blood, we used ReFACTor, a reference-free, unsupervised method that utilizes principal component analysis to estimate and adjust for cellular heterogeneity (Jaffe and Irizarry, 2014;Rahmani et al., 2016).In our study, we utilized the top five components identified by ReFACTor to control the potential effects of cell-type heterogeneity in EWAS analysis.

Gene expression data preparation
A subsample of 12 MZ twin pairs were included in the gene expression analysis.Briefly, total mRNA was extracted from whole peripheral blood.Then the RNA-Seq library was constructed and sequenced to obtain the sequenced data, which was mapped to the human genome by TopHat 2 (Kim et al., 2013).The gene expression level was evaluated by FPKM value through Cufflinks.

Epigenome-wide association analysis
For single CpG analysis, generalized estimating equation (GEE) model was applied to estimate the relationship between DNAm and HGS, with taking the correlation within each twin pair into account and adjusting for sex, age, and cell-type composition.We conducted this analysis using the geeglm function of geepack package in R software (version 4.1.0).The false discovery rate (FDR) was calculated to correct for multiple testing and genome-wide significance was defined as FDR <0.05.We annotated the identified genomic CpGs to the nearest genes by biomaRt package in R software.Causal inference.

Causal inference
The potential causal relationships between genome-wide significant CpGs and HGS were estimated using the ICE FALCON method, a regression-based method for causal inference in twin pair or family design (Li et al., 2005;Li et al., 2020).Estimations of β self , β co-twin and β' self , β' co-twin were calculated by the GEE model, where β co-twin is the estimation of family confounding proportion; β self represents the overall correlation including family confounding and casual proportion; and β'cotwin and β' self were the estimations of the full model.If the absolute difference between β self and β' self is greater than the absolute difference between β co-twin and β' co-twin (ratio >1.5), it suggests a causal relationship.However, if the association is a result of family confounding, the absolute difference between β self and β' self would be similar to the absolute difference between β co-twin and β' co-twin .

Differentially methylated region (DMR) analysis
Considering that DNAm of adjacent CpGs may be functionally or spatially linked, exploring genomic regions containing biologically linked DNAm may provide new perspectives compared to the analysis of single CpGs.Hence, we explored the DMRs potentially related to HGS by comb-p (Pedersen et al., 2012).The significantly enriched DMRs were identified using the Stouffer-Liptak-Kechris (slk) correction method with a significance threshold of p < 0.05.

Biological pathway analysis
We performed genomic region enrichment analyses using the Genomic Regions Enrichment of Annotations Tool (GREAT) (McLean et al., 2010).The CpGs that were identified with a p < 0.05 were submitted to the GREAT online platform for ontology enrichment.Annotation was conducted based on Genome Reference Consortium Human Build 37, employing the default "basal plus extension" association rule.Enrichment items with FDR<0.05 were deemed statistically significant.

EWAS power estimation
According to a recent computer simulation study on the power of EWAS based on twin designs, for traits with a heritability of 0.6 or higher, the maximal sample size of 63 pairs is required for statistical power exceeding 80% in utilizing trait-discordant twin design (Li et al., 2018;Wang et al., 2023).Compared to ordinary case-control design, twin design can significantly reduce the required sample size.Therefore, our sample size could meet the criteria for statistical power exceeding 80%.

DNA methylation and gene expression analysis
The correlation between top CpGs methylation level and corresponding gene expression level was calculated using Spearman's rank correlation.In addition, we also queried the position of significant CpGs in NCBI to illustrate the gene structures where these CpGs are located.

Quantitative methylation analysis of FBLN1 and ABHD14B
To validate CpGs located at FBLN1 and ABHD14B, we recruited 117 participants with low HGS (defined as male<28 kg, female<18 kg) and 117 healthy controls from the community.The method for HGS testing is consistent with the aforementioned approach.Participants attended interviews and venous blood samples were collected for DNAm level measurements.Primers for the FBLN1 and ABHD14B genes were designed to cover the majority of the CpGs identified in EWAS.The cleavage products were analyzed using MALDI-TOF mass spectrometry based on MassARRAY System (Bio Miao Biological Technology, Beijing, China).Subsequently, the resulting spectra were processed using the MassARRAY EpiTYPER software (Agena Bioscience, San Diego, California) to determine the methylation ratio.The Wilcoxon rank-sum test was used to compare the DNAm levels of GpGs between the two groups.The association between each CpG and low HGS was evaluated by logistic regression, with adjusting for age, sex, and BMI.The significance level was defined as p < 0.05.

Results
The epigenome-wide association study included 66 MZ with a median HGS of 32.0 kg (95% range 16.7, 57.0) and the median of absolute values of intra-pair HGS was 3.9 kg (95% range 0.3, 13.3).Significant correlations were observed between various clinical indicators within twin pairs, such as height, weight, BMI, fat rate, systolic and diastolic blood pressure, glucose, serum uric acid, cholesterol, triglycerides, HDL cholesterol, and LDL cholesterol, suggesting that utilizing a trait-discordant monozygotic twin design provided significant advantages (Supplementary Table S1).nine CpGs located at six different genes, including ERO1L, MFSD6, TMEM233, TYW1B, PAX8, and MRPL23.Methylation level of 16 CpGs (located at FBLN1, ABHD14B, MFSD6, TMEM233, and PAX8 genes) were positively associated with HGS, indicating that hypermethylation of these CpGs in twins with higher HGS.While nine CpGs (located at RXRA, ERO1L, TYW1B, and MRPL23 genes) were negatively associated with HGS, indicating that hypomethylation of these CpGs in twins with higher HGS (Table 1).

Causal inference analysis
Table 2 presents the estimation of the causal inference between genome-wide significant CpGs and HGS.In the causal inference from DNAm to HGS, the causal effect was supported for 16 CpGs located at FBLN1, RXRA, ABHD14B, MFSD6, and TYW1B.The causal effects from HGS to DNAm were also observed among two CpGs located at FBLN1 and RXRA.

Region-based analysis
A total of 21 DMRs were detected for HGS (Table 3).Notably, among these DMRs, 11 DMRs (A, C, E, F-J, L, S, and T) indicated a positive correlation with HGS, whereas five (B, D, Q, R, and U) near/ at CPLX1, PLEC, AKR1D1, WNK2, and ZNF597 exhibited a negative association.In the case of five DMRs (K, M-P), the methylation levels were inconclusive (Figure 2).

Biological pathway analysis
The GO-enriched terms mainly highlighted calcium channels, neuromuscular junctions, and skeletal muscle organ development.Moreover, the pathway terms associated with HGS were significantly enriched, including calcium transport, choline biosynthesis, genes involved in the integration of energy metabolism, PKC-catalyzed phosphorylation of inhibitory phosphoprotein of myosin phosphatase, and the Notch signaling pathway (Table 4).Some items are related to skeletal muscle function in previous studies and have shown a potential role of epigenetic changes in the denervation of calcium channels and neuromuscular junctions.

DNA methylation and gene expression analysis
A total of 12 twin pairs with a median age of 53 years (95% range 43-65) and a median HGS of 36.9 kg (95% range 18.4-61.3)were included in the analysis.Supplementary Table S2 presented the results of DNA methylation and gene expression analysis.
Methylation level of CpGs located at FBLN1, RXRA, and MRPL23 gene were positively associated with corresponding gene expression level.Methylation level of CpGs located at ABHD14B were negatively associated with gene expression.As shown in Supplementary Figure S1, all of these CpGs are located within gene bodies.

Quantitative methylation analysis of FBLN1 and ABHD14B
We quantified nine CpGs (FDR< 0.05) mapped to FBLN1 using the Sequenom MassARRAY platform.Six CpGs were negatively associated with HGS weakness (Supplementary Table S3).Four out   of the five CpGs mapped on ABHD14B were quantified and one CpG (chr3: 52,007,647) were not detected in the validation experiment.Just one CpG (chr3: 52,007,595) was validated to be hypomethylated in the low HGS group (Supplementary Table S3), although this CpG did not reach genome-wide significance in the EWAS analysis (FDR = 0.054).

Discussion
In the present study, we detected epigenetic variants of handgrip strength using EWAS based on monozygotic twins.Our analysis identified several CpGs, genes, DMRs, and pathways that might elucidate the mechanism of muscle strength change.As additional validation, several candidates CpGs mapped to FBLN1 and ABHD14B were quantified and validated.
Research has indicated that the effects of DNA methylation are associated with its occurrence at different positions in the genome and DNA methylation occurring in gene bodies may be related to gene expression (Aquino et al., 2018).Similarly, we found that the methylation levels of CpGs within FBLN1, RXRA, and MRPL23 are positively correlated with gene expression.The methylation of CpGs in FBLN1 is also positively associated with HGS.Fibulin-1, encoded by FBLN1, is an important protein in the structure of elastic fibers and basement membranes of various tissues.Previous studies suggested that Fibulin-1 implicated in tissue organogenesis during myotome development, bone formation, ossification, and digits of the development limbs (de Vega et al., 2009;Bohlega et al., 2014;Cooley et al., 2014).Methylation levels of CpGs located at RXRA and MRPL23 were negatively associated with HGS in EWAS.RXRA is a member of the nuclear receptor superfamily involved in lipid, glucose, energy, and hormone metabolism.Evidence from animal experiments indicates that RXRA promotes adipogenesis and accelerate fat accumulation (Pan et al., 2023).Previous research reported that the DNAm level at four CpGs located at RXRA (chr 9: 136,355,569-136,355,885) was negatively associated with bone mass, which might affect muscle strength (Harvey et al., 2014;Kim et al., 2018).The mechanism underlying the relationship between MRPL23 and muscle strength remains unclear at present and awaits further investigation in future studies.The methylation levels of CpGs located within the gene body of ABHD14B is positively associated with HGS and negatively correlated with gene expression.A previous study found that cortical ABHD14B levels were negatively associated with motor resilience, which reflects muscle endurance and motor function in the elderly (Buchman et al., 2023).
The present study also identified 21 DMRs in genomic regions located near or at 21 genes, of which FBLN1, SLC30A8, CST3, and  (Sprouse et al., 2014).The protein encoded by CST3 has been implicated as a biomarker for sarcopenia, as measured by the creatinine-to-cystatin C ratio (Tabara et al., 2020;Tabara et al., 2021).PLEC encodes a cytoskeletal protein that maintains tissue integrity by regulating intracellular signaling in response to mechanical stimulation (Rice et al., 2019).SOCS3 has been identified as a crucial mediator in myogenesis.SOCS3 overexpression during myogenesis enhanced myogenin and αactin mRNA expression (Caldow et al., 2011).The protein encoded by CTSD maintains cellular protein homeostasis and is regarded as a bone-related biomarker of osteoporosis (Deng et al., 2021;Wang et al., 2022;Wu et al., 2022).
Our findings provide evidence of causation underlying DNAm-HGS.We observed HGS might be a response to the methylation level of several CpGs located at different genes.The potential function of FBLN1, RXRA, and ABHD14B in muscle strength has been discussed above.Recent research has indicated that the methylation level of RXRA promoters could change in response to the external stimulus of mechanical loading (Krstic et al., 2022), which might explain the reason for the DNAm response to HGS changes.Additionally, studies have reported an association between  lower plasma vitamin D levels and elevated RXRA gene methylation (Curtis et al., 2019;Krstic et al., 2022).Considering the role of vitamin D in muscle strength, the changes in HGS-induced RXRA DNA methylation alterations may also be attributed to the potential confounding effect of vitamin D.
The strengths and advantages of the current study are as follows.First, our study included monozygotic twin pairs, which controlled the genetic background and enhanced the credibility of the results.Second, we investigated the causal relationship between DNAm and HGS, providing a deeper understanding of the effect of epigenetic modifications on HGS.Thirdly, our study is the first EWAS of HGS in the Asian population.Given the diverse genetic backgrounds and environmental exposures among different ethnic groups, our research illuminates the intrinsic physiological mechanisms contributing to the variation in HGS.
There are limitations of our study that must be acknowledged.First, our study did not consider physical activity and exercise due to a lack of relevant information, which can directly affect HGS or DNA methylation.The second limitation lies in the use of DNA methylation data derived from blood rather than the more intuitively relevant muscle tissue.Despite the challenges in obtaining muscle tissue, we acknowledge the tissue specificity of DNA methylation.Previous study has also indicated differences in the epigenetic feature between muscle tissue and blood (Slieker et al., 2013).However, a study conducted by Pilling et al. identified potential gene expression characteristics related to muscle strength in the blood and revealed several genes previously reported to be associated with muscle strength (Pilling et al., 2016).Additionally, some studies suggest that while the use of whole blood DNA methylation data may only reflect changes in blood, it can still provide information about variations in muscle mass (Livshits et al., 2016).Therefore, the use of blood samples for epigenetic studies on muscle strength may be a feasible approach (Soerensen et al., 2019).Finally, we recognize that HGS serves as an indicator of upper limb strength rather than a comprehensive measure of whole-body strength.Hence, future research should employ more representative indicators to validate our findings.

Conclusion
In summary, our study utilizing monozygotic twin pairs detected many DNAm variants that may be associated with handgrip strength, particularly the loci within FBLN1 and ABHD14B.Our findings provide new clues to the epigenetic modification underlying muscle strength.study design, data collection and analysis, decision to publish, or preparation of the manuscript.

FIGURE 1 Circular
FIGURE 1 Circular Manhattan plot for the epigenome-wide association study of HGS.The number of chromosomes and the -log 10 of p-values for statistical significance are shown.The dots represent the observed CpGs.

FIGURE 2
FIGURE 2 Differential methylation patterns for the identified DMRs.11 DMRs (A, C, E-J, L, S, and T) indicated a positive correlation with HGS, whereas five (B, D, Q, R, and U) exhibited a negative association The vertical axis shows the coefficient for the association of each CpG with HGS, and the horizontal axis shows the chromosome positions with the black points indicating each CpG.The black line presents the methylation pattern for each DMR.BP, base pair; chr, chromosome.

TABLE 1
The results of the epigenome-wide association study on HGS (FDR<0.05).
FDR: false discovery rate; NA, not available.

TABLE 2
Causal inference analysis between CpGs and handgrip strength.

TABLE 3
The results of annotation to significant differentially methylated regions (DMRs) (slk corrected p < 0.05).

TABLE 4
The top GREAT ontology enrichment for regions potentially related to handgrip strength by using a binomial test.

TABLE 4 (
Continued) The top GREAT ontology enrichment for regions potentially related to handgrip strength by using a binomial test.