Epigenome-wide exploratory study of monozygotic twins suggests differentially methylated regions to associate with hand grip strength

Hand grip strength is a measure of muscular strength and is used to study age-related loss of physical capacity. In order to explore the biological mechanisms that influence hand grip strength variation, an epigenome-wide association study (EWAS) of hand grip strength in 672 middle-aged and elderly monozygotic twins (age 55–90 years) was performed, using both individual and twin pair level analyses, the latter controlling the influence of genetic variation. Moreover, as measurements of hand grip strength performed over 8 years were available in the elderly twins (age 73–90 at intake), a longitudinal EWAS was conducted for this subsample. No genome-wide significant CpG sites or pathways were found, however two of the suggestive top CpG sites were mapped to the COL6A1 and CACNA1B genes, known to be related to muscular dysfunction. By investigating genomic regions using the comb-p algorithm, several differentially methylated regions in regulatory domains were identified as significantly associated to hand grip strength, and pathway analyses of these regions revealed significant pathways related to the immune system, autoimmune disorders, including diabetes type 1 and viral myocarditis, as well as negative regulation of cell differentiation. The genes contributing to the immunological pathways were HLA-B, HLA-C, HLA-DMA, HLA-DPB1, MYH10, ERAP1 and IRF8, while the genes implicated in the negative regulation of cell differentiation were IRF8, CEBPD, ID2 and BRCA1. In conclusion, this exploratory study suggests hand grip strength to associate with differentially methylated regions enriched in immunological and cell differentiation pathways, and hence merits further investigations. Electronic supplementary material The online version of this article (10.1007/s10522-019-09818-1) contains supplementary material, which is available to authorized users.


Introduction
Aging is accompanied by a general loss of physical capacity that is often associated with adverse agerelated outcomes and thus may have significant influences on the quality of life. For example, agerelated decline in skeletal muscle strength, which is an important component of sarcopenia and frailty, is associated with physical disability, morbidity and mortality (Beaudart et al. 2017).
Hand grip strength is a simple and easily administered measurement of muscle strength that correlates well with overall muscle strength (Visser et al. 2000;Wang et al. 2005;Bohannon et al. 2012). In accordance, the average hand grip strength peaks in early adulthood, followed by a period of continuance, and then starts declining from midlife and throughout life (Frederiksen et al. 2006;Dodds et al. 2014). Hand grip strength has been associated with a number of adverse health outcomes, such as depression (van Milligen et al. 2011;Fukumori et al. 2015), cardiometabolic risk markers and diseases (Mainous et al. 2015;Lawman et al. 2016;Lee et al. 2016), and multimorbidity (Cheung et al. 2013;Yorke et al. 2015), as well as with an increased overall mortality risk (Newman et al. 2006;Arvandi et al. 2016). In a study of almost 140,000 individuals from 17 different countries, Leong et al. (2015) found associations between hand grip strength and incident cardiovascular disease and further confirmed that hand grip strength is a strong predictor of all-cause and in particular cardiovascular mortality. Recently, a study of 502,293 participants from the UK Biobank confirms an association to all cause mortality and cause specific mortality from cardiovascular disease, yet also puts forward associations to cause specific mortality from respiratory disease, chronic obstructive pulmonary disease and cancers, as well as to incidences of cardiovascular disease, respiratory diseases, chronic obstructive pulmonary disease and cancers (Celis-Morales et al. 2018). In addition, recent literature reviews firmly validate the associations between hand grip strength and cognitive decline (Fritz et al. 2017;Rijk et al. 2016) and functional status such as activity of daily living (Rijk et al. 2016).
Previous studies have estimated that the heritability of hand grip strength is approximately 50-60% (Frederiksen et al. 2002). However, very few consistent genetic associations have been found, and large genome-wide association studies (GWAS) of hand grip strength are scarce compared to other traits, thus leaving the molecular background of hand grip strength largely unaccounted for (Garatachea and Lucía 2013;Chan et al. 2015). The very first hand grip strength GWAS included 2088 adults (aged 55-85 years) and revealed no genome-wide significant results (Chan et al. 2015). Later, a meta-analysis of 27,581 individuals over 65 years of age from 14 different cohorts discovered one genome-wide significant association, rs752045, in an intergenic region on chromosome 8, which was replicated in 6363 individuals from three additional cohorts (Matteini et al. 2016). rs752045 is located in a DNase hypersensitivity site known to affect a binding motif of the CCAAT/ enhancer-binding protein-b (CEBPB), which has been implicated in muscle cell differentiation and muscle repair (Matteini et al. 2016). Recently, an even larger multi-center GWAS of 195,180 individuals identified 16 loci including genes involved in processes such as structure and function of skeletal muscle fibers (ACTG1), neuronal maintenance and signal transduction (PEX14, TGFA, SYT1), and monogenic syndromes with involvement of psychomotor impairment (PEX14, LRPPRC and KANSL1) (Willems et al. 2017).
In addition to genetic studies, blood-based studies of molecular signatures may prove valuable for determining systemic factors and pathways involved in muscle strength maintenance and age-related decline. A gene expression analysis of blood messenger RNA from 698 people found the expression of the CEBPB gene to be significantly associated with hand grip strength as well as with physical performance (Harries et al. 2012). In a larger transcriptome-wide association study (TWAS), Pilling et al. (2016) investigated the association between hand grip strength and gene expression levels in a total of 7781 individuals aged 20-104 years from 4 independent cohorts. They identified 208 differentially expressed genes, half of which were novel in muscle-related literature, thus warranting future work on mechanisms underlying these associations (Pilling et al. 2016). Recently, Murabito et al. (2017) examined the association of 229 whole-blood microRNAs (miRNAs) with hand grip strength in a population of 5668 individuals aged 24-90 years. The study identified 93 significant miRNAs, for which the potential gene targets were involved in biological pathways of importance for muscle and aging (e.g. metabolism) .
To our knowledge, only one study has compared blood-based DNA methylation with hand grip strength. In 2015, Marioni et al. performed an epigenome-wide association study (EWAS) of hand grip strength in 1085 individuals from the Lothian Birth Cohort (mean age of 69.5 years) and found no CpG sites with a p value \ 10 -5 . However, looking at the related phenotype skeletal muscle mass, a bloodbased EWAS in 1550 middle-aged women, including 50 monozygotic twin pairs discordant for hand grip strength, identified a number of significantly associated differentially methylated genomic regions, of which some showed replication in 1196 individuals: these included the DNAH12, ZFP64, TP63 and GUSBP11 genes, which have previously been linked to muscle differentiation and function (Livshits et al. 2016). Hence, the latter study suggests the importance of epigenetic approaches for understanding the molecular basis of muscle phenotypes such as hand grip strength, and furthermore supports the use of blood as an informative tissue in studies of muscle function.
With the aim of exploring the association between hand grip strength and DNA methylation profiles in blood cells, we studied a sample of 672 monozygotic twins with an age span of 55-90 years. The study took advantage of the twin data by applying an intra-pair approach, controlling the influence of genetic variation and lowering the influence of potential confounders such as differences in rearing environment. Furthermore, we performed a longitudinal EWAS in 192 elderly individuals (age 73-90 years at intake), for whom data on hand grip strength were available from up to five survey waves conducted over 8 years.

Study populations
The study population comprised 672 monozygotic twins drawn from population-based, nation-wide surveys conducted within the framework of the Danish Twin Registry (Skytthe et al. 2013). Of these, 480 were participants from the Study of Middle Aged Danish Twins (MADT) and 192 were participants from the Longitudinal Study of Aging Danish Twins (LSADT) (Skytthe et al. 2013).
MADT was initiated in 1998 as a Danish nationwide study of 4314 twins randomly selected from birth cohorts spanning 1931-1952. In 2008-2011 a followup study was conducted of all eligible twin pairs originally enrolled (Skytthe et al. 2013). The present study included all intact monozygotic (MZ) twin pairs from the follow-up study. LSADT is a cohortsequential study of Danish twins aged 70 years or more, initiated in 1995. Surviving twins were followed-up every second year until 2005 (McGue and Christensen 2007). In 1997, whole blood samples were collected from 689 same-sex twins and the present study included all MZ pairs from this wave. For replication purpose we used data on 275 participants from the Study of Birth Weight Discordant Twins (BWD), a nation-wide study conducted in 2008-2010 of the most birth weight discordant twins in the twin registry (Frost et al. 2012). These samples were previously analyzed for epigenome-wide association with birth-weight discordance by Tan et al. (2014), who reported no significant findings. The demographics of the study sample are shown in Table 1.
Informed consents were obtained from all participants and all surveys were approved by The Regional Committees on Health Research Ethics for Southern Denmark (S-VF-19980072, S-VF-20040241 and S-20090033) and were conducted in accordance with the Helsinki II declaration.

Phenotype data
Hand grip strength was measured by a hand-held dynamometer (SMEDLEY'S dynamometer, Scandidact, Kvistgaard, Denmark), using a maximum of three measurements with the strongest hand. For LSADT the measurements were performed in 1999, while the blood sample was drawn in 1997, giving a mean age difference of 1.79 years (standard deviation (SD) = 0.28). In analysis of the LSADT data, the covariate values (including age) were from 1997. For MADT and BWD all data were collected in the same survey wave, yet not necessarily on the same date. On average, the difference between blood sampling and measurement of hand grip strength was less than 2 weeks.
For the 192 LSADT individuals longitudinal data were available. Survey waves were performed in 1999, 2001, 2003, 2005 and 2007. In summary, 39 individuals took part in the baseline wave (1999) only, while 46, 36 and 34 individuals took part in 2, 3 or 4 waves, respectively, and 37 individuals took part in all 5 waves. The individual-level hand grip strength values decreased in the study population as a whole (see Online Resource 1, Supplementary Fig. 1), with an average decline from intake to the last measurement of -0.49 kg/year (SD 1.42, range -5.9; 7.3).
DNA methylation data DNA was isolated from buffy coat by salt precipitation, using either a manual protocol (Miller et al. 1988) or a semi-automated protocol with the Autopure System (Qiagen, Hilden, Germany). Bisulfite conversion of 500 ng genomic DNA was performed using the EZ Methylation Gold kit (Zymo Research, Orange County, CA, USA) and DNA methylation was measured using the Infinium HumanMethylation450K BeadChip (Illumina, San Diego, CA, USA) at Leiden University Medical Center or GenomeScan B.V., As the BWD is made up by two age groups, they are listed separately here No. number of, SD standard deviation Leiden, The Netherlands, according to the manufacturer's protocols. The BeadChip images were scanned using the iScan system. The DNA methylation data were prepared on four different occasions giving four datasets; the BWD dataset (N = 296), the MADT dataset (N = 498) and two separate LSADT datasets (N = 224 and N = 86, respectively). Quality control and normalization were performed similarly to Tobi et al. 2015 and separately for each dataset. Data pre-processing was done using the free R packages MethylAid (van Iterson et al. 2014) and minfi (Aryee et al. 2014). Samples were excluded, first if less than 95% of the probes had a detection P-value \ 0.01, or second if the samples failed based on the internal quality control probes used by MethylAid (van Iterson et al. 2014). The sex of the individuals was confirmed by considering the values of the X chromosome probes by multidimensional scaling. The DNA methylation data were annotated using GRCh37/hg19. In total, 4 BWD twin pairs and 6 MADT twin pairs were excluded due to poor sample quality (8) or sex discrepancy (2). Probes were excluded if they had a detection P-value [ 0.01, a raw intensity value of zero, a low bead count (\ 3 beads), a measurement success rate below 95% or were identified as being cross-reactive (Chen et al. 2013). In total, between 32,055 and 33,316 CpGs out of the 485,512 CpGs on the array were excluded for of the four datasets. After probe filtering, a sample success rate of 0.95 was applied; no samples were excluded. Normalization was carried out by Functional normalization (Fortin et al. 2014) with four principle components (PCs) for the MADT and LSADT datasets and 5 PCs for the BWD dataset. The present study includes CpGs overlapping between the MADT and LSADT datasets (452,178 CpGs). Before statistical analysis, the b-values were transformed to M-values by logit transformation, presenting a distribution more suitable for statistical testing (Du et al. 2010).
Blood leukocyte subtypes (monocytes, lymphocytes, basophils, neutrophils and eosinophils) counted using a Coulter LH 750 Haematology Analyser were available for 483 of the 492 MADT individuals and for 212 of the 292 BWD individuals. For the remaining individuals the blood counts were imputed by a modified version of the PredictCellComposition method (www.github.com/mvaniterson/ predictcellcomposition), employing partial least squares regression after setting values of 0 to 0.0001, i.e. inflating percentages below the minimum value observed (0.0099), and transforming the cell composition using an isometric log ratio transformation.

Cross-sectional analyses
A principle component analysis (PCA) of the methylation data found the top 4 PCs to describe the majority of the variance in the data, as the subsequent PCs accounted for less than 3% of the variance (data not shown). Still, as PC1 correlated highly (correlation coefficient: 0.99) with sex, only PC2-4 were included in the statistical models. Furthermore, as the technical co-variates DNA plate number, 450K array number and dataset reflected much of the same variance, only the dataset co-variate was included. Cell counts for monocytes, lymphocytes and eosinophils were included in all statistical models, while basophil counts were excluded due to low variability, and neutrophil cell counts were excluded as neutrophil and lymphocyte cell counts reflected much of the same variance. Finally, age and height were included in all statistical models, as they associate with hand grip strength (Frederiksen et al. 2006), which was also found for the present data (data not shown).
Cross-sectional EWAS with single CpG methylation levels as the outcome and hand grip strength as the explanatory variable, were performed with the use of linear regression models, first on an individual level (individual level EWAS) and next on a twin pair level using an intra-pair design (twin pair level EWAS).
The individual level analysis was carried out with a mixed effect model using twin pairing as a random factor to accommodate for correlations between twins in a pair. Sex, age at blood sampling, height, cell counts (monocytes, lymphocytes and eosinophils) and PC2-4 were included as fixed effects and a dataset variable as a further random effect.
In the twin pair level analysis, a linear regression model was applied (similarly to Starnawska et al. 2017) to test the association between intra-pair DNA methylation difference and intra-pair difference in hand grip strength, while adjusting for sex, mean age of the twin pair at time of blood sampling, as well as for the intra-pair difference in cell counts (monocytes, lymphocytes and eosinophils), height and PCs 2-4.
For the 321 intact twin pairs of the present study population the average hand grip strength was 34.0 kg (SD 12.2, range 7-70), while the average intra-pair difference was 4.4 kg (SD 3.8, range 0-21). The intrapair difference for a given co-variate was calculated by subtracting the co-variate value of the twin with the lower hand grip strength from the co-variate value of the co-twin with the higher hand grip strength. In a sub-analysis, the twin pair level model was applied to the 50% most discordant twin pairs with respect to hand grip strength, thus removing twin pairs with no or little intra-pair difference.
The replication analysis in the 275 BWD individuals was performed as described above, with additional adjustment for birth weight or intra-pair difference in birth weight.
All statistical analyses were carried out in R version 3.3.1.

Longitudinal analyses
An EWAS of hand grip strength over time was carried out by a growth curve model, which employs a linear mixed-effects model with hand grip strength at the different time points as the outcome and the single CpG methylation level, time from baseline to survey, sex, age at blood sampling, cell counts (monocytes, lymphocytes and eosinophils), height and PCs 2-4 as fixed effects. A dataset variable, a twin pairing variable and a twin ID variable (allows individual varying baseline levels) were included as random effects, whereby the twin ID specific effect was nested within the twin pair specific effect. This initial random intercept model was extended by a random slope model (allows individual varying changes in hand grip strength over time). Considering that the two models gave very similar results, although slightly higher p-values in the random slope model (data not shown), only the results from the random intercept model are reported. All statistical analyses were carried out in R version 3.3.1.

Pathways analyses of EWAS findings
CpGs with an uncorrected p-value \ 10 -3 were annotated using the information by Illumina and the corresponding genes were used for overrepresentation enrichment analysis (ORA) against the genes on the 450K array using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome databases available through the WEB-based GEne SeT AnaLysis Toolkit (WebGestalt) (Wang et al. 2017). Moreover, an ORA sub-analysis restricted to the 134,070 CpG sites annotated to the Function Feature Type TSS200 (transcription start site (TSS) to -200 nt upstream of TSS) or TSS1500 (-200 to -1500 nt upstream of TSS) was performed. The number of genes investigated and annotated to the databases can be seen in the Online Resource 1, Supplementary Table 1.

Differentially methylated regions
In order to identify differentially methylated regions (DMR), the comb-p approach was applied (Pedersen et al. 2012) for the results obtained in the crosssectional and longitudinal EWAS. The comb-p has been reported to have high control of false-positive rate and the best sensitivity compared to other popular DMR tools, such as bumphunter, DMRcate and probe lasso (Mallik et al. 2018). The following settings were used: p-value was set to 0.1 to seed a region; maximum 500 bps to skip before finding next seed, if new seed was not found, the region is ended. Gene regions identified were annotated using the BioMart (www. biomart.org) and Ensemble (www.ensembl.org) databases. The gene regions passing correction for multiple testing by false discovery rate (FDR) (see below) were investigated by MSigDB genomic regions enrichment analysis (GREA) by the Genomic Regions Enrichment of Annotations Tools (GREAT) (http:// great.stanford.edu) using the BED file output from comb-p as input. The hypergeometric test over genes was applied.

Adjustment for multiple testing
In all analyses adjustment for multiple testing was performed by the Benjamini-Hochberg FDR correction method (Benjamini and Hochberg 1995). The suggested genome-wide significance threshold for EWAS analyses using 450K data has been proposed to be p \ 10 -6 (Tsai and Bell 2015). However, as we consider this an exploratory study, and in order to give an overview of the top results, we report all findings with an unadjusted p \ 10 -5 for the EWAS results.

Results
The demographics of the study sample are shown in Table 1. All Supplementary results can be found in the Online Resource 1.
Single CpG epigenome-wide association studies of hand grip strength Epigenome-wide association studies (EWAS) were performed both on the intake hand grip values for the entire study population, as well as on longitudinal hand grip measurements, which were available for the elderly twin sample ( No genome-wide significant associations were observed in the cross-sectional analyses of the 672 individuals (see Table 2A-C). Five CpG sites displayed a p-value below the cut-off of 1 9 10 -5 in the individual level EWAS: four of these were located in gene regions, i.e. in the D-glutamate cyclase (DGLUCY), NLR family CARD domain containing (NLRC5), tec protein tyrosine kinase (TEC) and iroquois homeobox 4 (IRX4) genes, respectively, and 1 CpG site was in an intragenic region, the nearest gene being MT-RNR2 Like 1 (MTRNR2L1). Similarly, the twin pair level EWAS displayed 5 CpG sites with a p-value \ 1 9 10 -5 : four were located in gene regions, i.e. in the solute carrier organic anion transporter family member 3A1 (SLCO3A1), calcium voltage-gated channel subunit alpha1 (CACNA1B), RNA binding motif protein 33 (RBM33) and LIM domain binding 2 (LDB2) genes, respectively and 1 CpG site was in an intragenic region, the nearest gene being Lengsin (LGSN). The CpG site cg01763947 in SLCO3A1 was also found to be associated in the analysis of the 50% most discordant twin pairs, together with 6 additional CpG sites.
None of the 16 top-ranked CpG sites listed in Table 2A-C replicated convincingly in the BWD twins (see the Online Resource 1, Supplementary Table 6).
In the overrepresentation enrichment analyses (ORA) against the KEGG and Reactome databases of the CpG sites with an uncorrected p-value below 1 9 10 -3 , none of the identified pathways reached statistical significance after correction for multiple testing. All ORA results with an uncorrected p-value \ 0.05 are summarized in the Online Resource 1, Supplementary Table 7.

Analysis of differentially methylated regions
The EWAS results obtained were used for identification of differentially methylation regions (DMRs) by the comb-p algorithm. In total, 70, 58, 51 and 63 DMRs displayed a FDR p-value below 0.05 for the individual level EWAS, the twin pair level EWAS, the twin pair level EWAS for the 50% most discordant pairs and the longitudinal EWAS, respectively (see Tables 3, 4, 5 and 6). Of these, 60, 49, 41 and 54 DMRs, respectively, were annotated to gene regions, while the remaining DMRs were annotated to intergenic regions.
Applying the comb-p output for functional enrichment analysis using the GREAT tool (http://great. stanford.edu), revealed 8 significant pathways for the individual level EWAS; of these pathways, 7 were related to immune function or (auto)immune diseases, while 1 was related to C-MYC transcriptional repression (see Table 7). No pathways passed correction for multiple testing for the twin pair EWAS, the twin pair level EWAS for the 50% most discordant pairs and the longitudinal EWAS. All GREAT results with an uncorrected p-value \ 0.05 are summarized in the Online Resource 1, Supplementary Table 8.

Discussion
In this study we explored the association between variation in the epigenome and physical functioning represented by hand grip strength using both crosssectional and longitudinal measures of hand grip strength. The use of twins enabled us to explore epigenetic association independent of the genetic background of the individuals.    Although we did not identify any epigenome-wide significant associations of single CpG sites, some biologically plausible genes were observed (Table 2): the top CpG site of the longitudinal EWAS was located in COL6A1 encoding the alpha 1 subunit of type VI collagen. Collagens are extracellular matrix proteins, maintaining the integrity of various tissues. Mutations in COL6A1 have been reported to cause, among others, Bethlem myopathy (omim.org/entry/ 158810) and Ullrich Congenital Muscular Dystrophy 1 (omim.org/entry/254090) characterized by increased skeletal muscle weakness and joint stiffness in infancy. Furthermore, a top CpGs of the twin pair EWAS was located in the calcium voltage-gated channel subunit alpha1B (CACNA1B) gene, encoding the alpha-1B pore-forming subunit of an N-type voltage-dependent calcium channel. Voltage-dependent calcium channels are involved in calciumdependent processes, including muscle contraction, hormone or neurotransmitter release, gene expression and cell death. Mutations in CACNA1B have been reported for, among others, dystonia 23, characterized by involuntary muscle contraction (omim.org/entry/ 614860). Both genes have been reported to be expressed across a broad range of tissues, including whole blood and skeletal muscle (www.genecards. org).
Only one study on the epigenetics of hand grip strength has been published so far (Marioni et al. 2015) with no significant CpG sites or pathways reported.    ) have explored hand grip strength. None of these studies report overlapping genes among their top findings, and none of the suggestive genes of the present study were reported in these studies, neither among the top findings nor among the nominal significant genes reported by Harries et al. (2012), Chan et al. (2015), and Matteini et al. (2016). Furthermore, no significant pathways were identified in the ORA of the present study; 4 pathways did relate to muscle contraction, yet to cardiac and vascular smooth (blood vessel) muscle contraction, not skeletal muscle function (Online Resource 1, Supplementary Table 7). Of the previous studies, a TWAS (Pilling et al. 2016), a genome-wide miRNA study ) and a GWAS (Willems et al. 2017) reported pathway analyses. Due to differences in methodology a complete comparison was not possible, nevertheless, if grouping pathways by the overall hierarchical level of functional events as defined by the databases, some overlap between the present ORA and the previous studies was seen: pathways related to gene expression (Pilling et al. 2016), signal transduction, the neuronal system and cell cycle/death ), the immune system (Pilling et al. 2016 andMurabito et al. 2017), DNA repair (Willems et al. 2017), RNA metabolism ) and protein metabolism (Pilling et al. 2016and Willems et al. 2017. Such overlap could indicate that these processes might relate to hand grip strength, although it must be stressed that these overall categories are broad and hence can reflect a rather broad array of biological functions. Contrary to the single CpG analysis, the analysis of differentially methylated regions (DMRs), with enriched power by combining multiple CpGs, did show significant findings; of the 203 unique gene regions identified (Tables 3, 4, 5, 6), 5 have been reported by previous studies on hand grip strength: the PLEKHB1 pleckstrin homology domain containing B1 (PLEKHB1) and O-6-methylguanine-DNA methyltransferase (MGMT) genes in the GWAS of Matteini et al. (2016) and Willems et al. (2017), respectively, and the tachykinin receptor 1 (TACR1), G protein-coupled receptor 19 (GPR19) and prostaglandin F2 receptor inhibitor (PTGFRN) genes in the TWAS by Harries et al. (2012). Interestingly, mutations in PLEKHB1 have been reported to cause amyotrophic lateral sclerosis, characterized by problems with muscle control and movement (omim.org/ entry/105400).
Moreover, a number of genes were observed in more than one comb-p analysis (Online Resource 1, Supplementary Table 8): the major histocompatibility complex, class II, DP alpha 1 (HLA-DPA1),  Major histocompatibility complex genes were also revealed in the gene sets passing correction for multiple testing in the genomic regions enrichment analysis (GREA) ( Table 7): the major histocompatibility complex, class I, C (HLA-C), major histocompatibility complex, class I, B (HLA-B), major histocompatibility complex, class II, DM alpha (HLA-DMA) and major histocompatibility complex, class II, DP beta 1 (HLA-DPB1) genes did, together with the (non-muscular) myosin heavy chain 10 (MYH10), endoplasmic reticulum aminopeptidase 1 (ERAP1), and interferon regulatory factor 8 (IRF8) genes, constitute the 7 pathways related to the immune system or (auto)immune diseases. Some of these genes have been reported to be mutated in diseases related to the immune system, e.g. HLA-DMA: rheumatoid arthritis and lupus, HLA-C: psoriasis and HIV, HLA-B: ankylosing spondylitis (arthritis of the spine) and IRF8 and HLA-DPB1: immune deficiencies (www. genecards.org). Moreover, when considering all nominally significant findings from the GREA (Online Resource 1, Supplementary Table 8), 31 out of 80 unique gene sets were related to immune function. These findings point to immune function as important for hand grip strength, which might not be surprising considering the well-known role of the immune system in muscle growth and regeneration (Tidball 2017), a process which declines during aging, contributing to the occurrence of sarcopenia (Domingues-Faria et al. 2016).
The last gene set found in the GREA related to negative regulation of cell differentiation and was composed of the IRF8, the CCAAT enhancer binding protein delta (CEBPD), inhibitor of DNA binding 2 (ID2) and BRCA1 DNA repair associated (BRCA1) genes. The MYC10, ID2 and BRCA1 genes are involved in cancer, while CEBPD is involved in juvenile diseases of the central nervous system, as well as in regulation of genes involved in immune and inflammatory responses (www.genecards.org). Furthermore, as mentioned in the introduction, an interaction partner of CEBPD, CEBPB, has been implicated in muscle cell differentiation and muscle repair (Matteini et al. 2016), and was reported in the TWAS by Harries et al. (2012). Also, ID2 has been related to muscle cell differentiation (www.genecards. org), indicating a role of these two genes in muscle biology, potentially mediated via repression of cell differentiation or immune function.  Finally, in the GREA the individual level EWAS displayed significant findings (Table 7), while the twin pair level EWAS and the longitudinal EWAS did not. One reason could be that the latter holds lower power, as it only included the elderly part of the study population. Moreover, the twin pair EWAS controls for genetic variation by analyzing intra-pair differences by linear regression, while the individual level EWAS accommodate the intra-pair correlation by inclusion of twin pairing as a random factor, which relies on the fulfillment of the assumptions of the model. Hence, it is possible that the individual level EWAS results are, least in part, influenced by genetic variations left unaccounted for.
A possible limitation of the present study is the use of DNA methylation data from whole blood rather than a tissue more intuitively relevant for a muscle related trait. It is well known that epigenetics is important for differentiation of and aging in muscle cells (Carrió and Suelves 2015;Zykovich et al. 2014), and that, although there is some overlap, differences in the epigenetic signature of muscle tissue and blood cells have been reported (Slieker et al. 2013). The findings by Livshits et al. (2016) on muscle mass do, however, indicate that the use of blood samples for epigenetic studies is a feasible approach for identification of muscle-related epigenetic changes, although obviously restricted to changes that are mirrored in blood. This line of thought is supported by the bloodbased gene expression study performed by Pilling et al. (2016) on hand grip strength, which found several genes previously reported to be associated to muscle function (Pilling et al. 2016). In any case, to decipher the muscle cell specific epigenetic signatures related to muscle strength, future studies on the less accessible muscle cell derived DNA methylation data and hand grip strength are highly relevant. Furthermore, as the age-related decline in muscle strength is a multifactorial and complex process affected by for instance cognitive, neural and hormonal processes (Manini and Clark 2012;Tieland et al. 2018), we cannot exclude that other physiological aging-related phenomena are reflected in the blood samples under study here, and that the results obtained might reflect  the test set such phenomena. Specifically, we can not exclude that the findings might, at least in part, reflect frailty, especially in the elderly part of the study population, as decreased hand grip strength is one component of the frailty phenotype (e.g. Li et al. 2019). Most importantly, this point is supported by the fact that the most significant biological pathways found in this study are highly implicated in immune function. In the literature, the strong connection between immune function and physical frailty has been extensively studied and discussed (Lu et al. 2016;Fuentes et al. 2017). In this sense, findings from this exploratory study might provide potential references for exploring the epigenetics of age-related loss of physical capacity.
In conclusion, in this exploratory study we identified several epigenetic markers of potential importance to physical functioning, as reflected by hand grip strength, using both cross-sectional and longitudinal data from twins. Furthermore, investigation of differentially methylated regions identified significant genomic regions and pathways related to the immune system and autoimmune diseases. However, additional studies are needed to further clarify the potential importance of these markers and regions. Such clarification could elucidate the mechanisms underlining maintenance and loss of muscle strength and age-related loss of physical capacity.