Contribution of Heritability and Epigenetic Factors to Skeletal Muscle Mass Variation in United Kingdom Twins

Context: Skeletal muscle mass (SMM) is one of the major components of human body composition, with deviations from normal values often leading to sarcopenia. Objective: Our major aim was to conduct a genome-wide DNA methylation study in an attempt to identify potential genomic regions associated with SMM. Design: This was a mixed cross-sectional and longitudinal study. Setting: Community-based study. Participants: A total of 1550 middle-aged United Kingdom twins (monozygotic [MZ] and dizygotic [DZ]), 297 of which were repeatedly measured participated in the study. MainOutcomeMeasure: Appendicularleanmassassessedusingdual-energyX-rayabsorptiometry technology, and methylated DNA immunoprecipitation sequencing DNA methylation profiling genome-wide were obtained from each individual. Results: Heritability estimate of SMM, with simultaneous adjustment for covariates obtained using variance

L ean body mass (LBM), in particular skeletal muscle mass (SMM) is one of the 3 major components of body composition, which also includes fat body mass and bone mass. As for the other 2 components, SMM is highly important for normal physiology and metabolism, and deviations from normal values are often associated with various pathological conditions, notably sarcopenia, particularly in women (1). Sarcopenia is as a rule defined as age-related reduction in muscle mass and muscle strength, and affects women more than men with the prevalence as high as 30% for those above 60 years old (2,3). In the elderly, the loss of LBM is correlated with profound physical impairment and disability with severe clinical consequences, including mobility loss, osteoporosis, increased fracture risk, dyslipidemia, insulin resistance, and increased mortality (4).
Despite its clinical significance and LBM has a strong familial component (5)(6)(7)(8), replicated and confirmed specific genetic polymorphisms are rare and can explain only a minor part of the muscular mass variation (9 -12). To date, no confirmed genotypes are associated with accelerated sarcopenia or frailty, or can explain a significant part of the interindividual variation in SMM.
Twin studies have made a remarkable contribution to the understanding of the molecular basis of human complex traits via modern high-throughput genetic and genomic analyses, and recently epigenetic approaches (13,14). There have been major advances in epigenetics related to carcinogenesis (15,16), psychiatric diseases (17), rheumatoid arthritis (18,19), and others. However, we are not aware of studies exploring the epigenetics of muscular mass variation, despite suggestions of its importance (20), and recent work has explored age-related differential methylation signatures in human skeletal muscle (21).
The main aim of this study was therefore to conduct a genome-wide DNA methylation study of a large sample of United Kingdom twins, using methylated DNA immunoprecipitation sequencing (MeDIP-seq) in an attempt to identify potential genomic regions where DNA methylation levels are associated with SMM variation, with a secondary aim to obtain good estimates of heritability.

Study sample
The data examined in the present study were from the TwinsUK Adult Twin Registry, described in detail elsewhere (22). The sample had been collected from the general population through national media campaigns in the United Kingdom without ascertainment for any of individual characteristics, diseases or traits. In the present study, 1196 individuals (119 dizygotic [DZ] and 428 monozygotic [MZ] twin pairs and 102 singletons, 56 MZ and 46 DZ twins without sibling measurements) were included having all measurements of interest. All studied individuals were females with age range between 17 and 82 years and average 51.8 Ϯ 13.7 years. The mean body mass index was 25.2 Ϯ 4.7 kg/m 2 , and LBM ranged between 26.0 and 60.9 kg, with mean LBM ϭ 39.9 Ϯ 5.4 kg. All participants gave written informed consent before entering the study and the St Thomas' Hospital research ethics committee had approved the project.

Muscle mass phenotype
All 3 major body composition components, ie, bone mineral density, fat body mass, and LBM, were measured by using standard whole body dual-energy X-ray absorptiometry method (5), following manufacturer's recommendations (QDR 4500W system; Hologic, Inc). Briefly, at installation, the manufacturer's engineer calibrated the instrument, and then daily quality control (QC) scans were performed using the spine phantom. Intrascanner reproducibility, expressed as a coefficient of variation from duplicate measurements in healthy volunteers 1 week apart, was 0.8% at the lumbar spine. Both twins within a pair were scanned on the same day. For the purposes of this study, we defined SMM as the sum of LBM measurements at the 4 limbs (appendicular lean mass), and not total LBM, which is biased by measurement of nonmuscular soft tissue, in particular viscera.

Smoking scores
The present sample included 1100 individuals for whom the information on smoking habits was available. Of these, in 446 individuals were smokers, which included 320 current smokers and 126 exsmokers.
For blood cell composition, whole-blood cell (WBC) subtype counts were obtained for 441 individuals from the replication sample (n ϭ 1196) using fluorescence-activated cell sorting (FACS) of peripheral blood (23). WBC subtype cell counts were available for 4 cell types: neutrophils, eosinophils, monocytes, and lymphocytes.

DNA extraction and MeDIP-seq
Whole-blood samples (6 mL) were collected and stored at Ϫ80°C in EDTA tubes. Genomic DNA was extracted using the BACC3 Genomic DNA Extraction kit (Nucleon) and stored in TE buffer at Ϫ20°C. A total of 1.5 g of genomic DNA was fragmented to a smear of 200 -500 bp with the Bioruptor NGS System (Diagenode) and was subsequently end repaired, adenylated, and adapter ligated using the Paired-End DNA Sample Prep kit (Illumina). Methylated DNA was immunoprecipitation using the Magnetic Methylated DNA Immunoprecipitation kit (Diagenode) as previously described (24). After efficiency and sensitivity assessment by quantitative PCR, MeDIP-seq libraries were prepared by amplification (Platinum pfx DNA Polymerase kit; Invitrogen), purification (Agencourt Ampure Beads; Beckman Coulter) and validation (Agilent BioAnalyzer analysis) followed by high-throughput sequencing (Illumina HiSeq2000) that generated approximately 50 million, 50-bp single-end reads.

MeDIP-seq DNA methylation quantification
After adapter and base quality trimming sequencing, reads were mapped to hg19 using BWA v0.5.9 (25). Alignments with low quality scores (Q Ͻ 10) and duplicates were filtered, which resulted in an average of 15 684 723 uniquely mapped reads that were subsequently extended to 350 bp to represent the average

Design of the study and statistical analysis
The analysis of the data included the following stages and corresponding methods and was based on 3 nonoverlapping subsamples, selected from the total available sample.
The first sample was used for identification of longitudinally stable DNA methylation regions for downstream association analysis. For the selection of these methylation regions, we computed Pearson correlations between the all-available longitudinal measurements of the methylation within the individual. The bins displaying significant intraindividual correlations were considered as longitudinally stable, and from now on will be called longitudinally stable bin (lsBIN). We assumed that the lsBINs would also show significant correlation between twins, which could be regarded as an additional confirmation of their nonsporadic nature. To this end, 292 individuals having sequential measurements of methylation with the same processing batch were selected from the entire available sample (385 individuals having at least 2 methylation measurements each). Five of them had 2 pairs of appropriate sequential measurements each (4 individuals had 4 measurements and one 3 measurements with the same batch), the rest 2 measurements only. The measurements were taken with at least 3 years' time interval (mean, 7.0; SD, 1.2; range, 3.4 -10.8). For the intrapair correlation, 104 pairs of MZ twins (219 pair of measurements) and 63 pairs of DZ twins (129 pair of measurements) with methylation measured at the same date and within the same batch were selected from the same total sample.
The second sample was used for initial identification of the bins potentially associated with SMM variation. For this aim, we selected 50 of the most discordant SMM pairs of MZ twins, implementing a simple formula dSMM i ϭ (dSMM i1 Ϫ dSMM i1 )/ 0.5(dSMM i1 ϩ dSMM i1 ). The 50 dSMM i were explored for the corresponding methylation differences, dMTL i . Only lsBINs were used in this analysis. The DXA scans and methylation were taken with less than or equal to 12 months between both types of assessment.
The third sample was used for validation in an extended dataset. At this stage, lsBINs identified above were examined in the remaining sample including 1196 individuals in total (see above). We tested correlation between SMM and methylation level and used multiple regression analysis where SMM was considered a dependent variable, methylation level, age, and smoking as covariates. In a subset of this sample (n ϭ 441) with available WBC, we also considered WBC subtype proportions for association with DNA methylation levels at the 7 top-ranked SMM-associated methylation signals and found no differences in the pattern of their association.
To evaluate the contribution of the additive genetic factors (heritability estimates, h 2 ) to variation and covariation of SMM and methylation levels of the significantly associated bins, we carried out uni-and bivariate variance decomposition analysis, based on a classical polygenic concept of the quantitative trait inheritance. The analysis was conducted by MAN statistical package for familybased samples (http://www.tau.ac.il/ϳidak/hid_MAN.htm).

Functional genomic and CpG island (CGI) annotations
Annotations were obtained from University of California Santa Cruz (UCSC) (http://hgdownload.soe.ucsc.edu/goldenPath/hg19/ database/cpgIslandExt.txt.gz). CGI shores were subsequently defined as 2-kb regions up-and downstream of CGI boundaries. Promoter (active, weak, and poised) and enhancer (strong and weak) chromatin states predicted by ChromHMM (27) for GM12878 (a lymphoblastoid cell line derived from blood) (28) were also obtained from UCSC (http://hgdownload.cse.ucsc. edu/goldenPath/hg19/encodeDCC/wgEncodeBroadHmm/ wgEncodeBroadHmmGm12878HMM.bed.gz). Distribution of the methylation signals by the levels of their longitudinal stability, and extent of twin correlations among CGIs, promoters, and enhancers was compared using standard ANOVA and t tests.

Heritability of SMM
The crude SMM measurements showed modest but statistically significant inverse correlation with age (r ϭ Ϫ0.097, P ϭ .006). The intraclass correlations of MZ and DZ twins for the age-adjusted SMM were high and significant: R MZ ϭ 0.799, P ϭ .0001 and R DZ ϭ 0.366, P ϭ .0008, respectively, suggestive of strong genetic influence. Indeed the heritability estimate obtained using variance decomposition analysis yielded h 2 ϭ 0.809 Ϯ 0.050.

Identification of longitudinally stable methylation signals
First, from the total 11 524 145 bins quantified genome-wide, those displaying zero methylation levels in more than 20% of the individuals were excluded, leaving 6 501 931 bins (56.4%) for further analysis (Table 1). Of these, only a minor portion of 723 029 bins (6.3% of the initial 11 524 145 bines) showed significant positive correlation between the longitudinal MeDIP-seq measurements within individuals, ranging between 0.114 and 0.905, with nominal P Ͻ .05.
Next, we computed intrapair correlations between all 6 501 931 bins, for MZ and DZ pairs separately, ie, R MZ and R DZ . To explore whether the correlation between the twins depends on the bin longitudinal stability, we computed the correlation through all selected bins on each chromosome between the repeated measurements (R repeated ) and the correlation coefficients between the twin pairs for the corresponding bins, according to their zygosity (R MZ and R DZ ), ie, R repeated was contrasted with R MZ (or R DZ ).
The correlation between the R MZ and R repeated was at a genome-wide average 0.51 (ranging between 0.41 and 0.58 per chromosome) and was consistently greater than the corresponding correlation between the R DZ and R repeated , which was at a genome-wide average of 0.38 (and varied from 0.32 to 0.41 per chromosome). The corresponding results for each chromosome are provided in Table 1, and exemplifying scatterplots for chromosomes of different size are shown in Figure 1, A and B. Figure 1C demonstrates clear significant positive correlations between the R MZ and R DZ for the corresponding bins, with consistent tendency, R MZ Ͼ R DZ . This relationship, expressed as R MZ /R DZ (considering bins with positive and significant R DZ ), shows substantial correlation (0.46 -0.53, depending on chromosome) and highly significant (P Ͻ .0001) correlation with R repeated (Figure 2). This suggests that genomic regions with evidence for genetic heritability are more likely to be longitudinally stable.
The statistical significance (P value) of the correlations of R repeated with R MZ (as well as with R DZ ) per specific bin varied widely, depending on bin and chromosome. However, when only bins with R repeated Ͼ 0 and P Ͻ .05 were selected, all the aforementioned correlations (ie, correlations between R repeated and R MZ , or R DZ ) became highly significant (P Ͻ 10 Ϫ8 ) for all chromosomes.

Identification of the methylation association with SMM variation in MZ twins
We tested whether the lsBINs were associated with SMM variation in the sample of 50 SMM discordant MZ twins using paired t test to compare methylation levels per lsBIN. Of 723 029 bin comparisons, 36 081 showed nominally significant (P Ͻ .05) association with SMM, which corresponds closely to 5% expected from type 1 error. These results could not be considered as entirely random, in particular for estimating the possible rate of false positive association signals. The present analysis was conducted on the small fraction of bins (6.3%), selected by QC and longitudinal stability tests. We therefore opted not to apply statistical correction technique to estimate the rate of false discovery, at this stage, as this would also increase the probability of type-2 error. Instead, of the above lsBINs, we selected only those that showed in paired t test for discordant MZ twins P Ͻ .01 and R repeated Ͼ 0.4. In total, 134 autosomal bins fitted the selection criteria, and were tested in our replication sample. The distribution of the P values per chromosome is shown on Supplemental Figure 1, and the lsBINs from the most significant t tests (P Ͻ .001) with their genomic positions are listed in Supplemental Table 1.

Replication of the association results in the extended population sample
At this stage, correlations of the methylation levels of the 134 aforementioned autosomal lsBINs with SMM (dependent variable) were examined in the remaining 1196 individuals. Results for 20 of them are presented on left hand side of Supplemental Table 2 as unadjusted Pearson  Table 2. We next tested whether the inclusion of the additional potential covariates, specifically age and smoking status (available for 1100 individuals) affected the aforementioned association results. The results are given on the right hand side of Table 2. Their inclusion did not change significantly the results, although age was a consistently nominally significant covariate for all the tested bins, whereas smoking did not show significant effects.
Whole blood is heterogeneous collection of cells and DNA methylation levels may reflect cellular composition. We therefore considered if the 7 top-ranked replicating SMM-association signals ( Table 2) were also associated with blood cell subtypes in a subset of individuals from the replication sample, where WBC data were available. We observed no evidence for association at nominal significance (P ϭ .05) between DNA methylation levels at these regions and the proportion of lymphocytes, neutrophils, basophils, and eosinophils.
We explored whether these SMM methylation results were likely caused by common genetic and/or environmental factors in SMM. We therefore conducted classical bivariate variance component analysis via estimation of contribution of pleiotropic genetic factors (genetic correlation, r G ) and shared environmental effects (environmental correlation, r E ) and summarized in Supplemental Table 3. It provides univariate estimates of heritability (contribution of additive genetic factors) to SMM and each of the 7 top methylation signals, and pairwise estimates of r G and r E between each of the bins and SMM. The heritability estimates for all bins were significant and varied between 0.271 Ϯ 0.181 and 0.753 Ϯ 0.091. Interesting, by Likelihood Ratio Test, genetic correlations for 5 pairs of 7 were statistically significant, with P values from 0.0477 to 0.0002. The remaining 2 bins on chromosome 21 (15253001-15253500) and 3 (189366751-189367250), showed nonsignificant r G and r E, when each of them were tested separately. However, by LRT, both correlations could not be constrained to zero (hypothesis of no correla- tion is rejected), with P ϭ .0068 and 0.0293, suggesting that statistical power of the sample is probably insufficient to choose between r G and r E . Other findings of interest concerned the environmental correlations (r E ) for DNA methylation, which were all statistically nonsignificant (P values ranged between 0.752 and 0.068), despite the fact that initially those bins were selected in the sample of discordant MZ twins.

Distribution of methylation signals among the CpG rich and poor genomic regions
Selected lsBINs (723 029 in total) relatively to initial total number of bins (11 524 145) were sorted as belonging to 1) CGIs, 2) CpG shores, or 3) not belonging to any of these 2 regions (CpG sea). Table 3 presents the summary of these data suggesting that there is a significant loss of lsBIN portion in CGIs (0.0017 vs 0.0032) vs clear enrichment in CpG shores (0.0525 vs 0.0285).
However, when the distributions of R repeated in these 3 regions were compared, we found that on average significantly (P ϭ 2.2 Ϫ21 -7.6 Ϫ16 ) higher values were observed in islands (0.1857 Ϯ 2.4E-03) in comparison with both CGI shores (0.1630 Ϯ 3.1E-04) and open seas (0.1666 Ϯ 8.0E-05). With respect to R MZ , there were no statistically significant differences between the CGIs and shores, but both had significantly (P ϭ 3.5 Ϫ36 -1.0 Ϫ300 ) higher R MZ in comparison with CpG seas.

Distribution of methylation signals among functionally different genomic regions: predicted enhancers and promoters
We conducted a comparative analysis of the distribution of lsBINs mapped to predicted enhancers and promoters (Table 4, see also Supplemental Table 4). We observed statistically significantly lower levels of R repeated in enhancers (average R repeated ϭ 0.1634 Ϯ 0.0004, n ϭ 53357) in comparison with promoters (average R repeated ϭ 0.1698 Ϯ 0.0013, n ϭ 8412) (t ϭ Ϫ7.84, P ϭ 4.7 Ϫ15 ), as expected. Within enhancers, the percentage of lsBIN decreases (from 8.56% to 2.29%) almost monotonically (Spearman's ϭ Ϫ1.0, P Ͻ .0001) with the average longitudinal stability (R repeated ) of the bin, that is, DNA methylation bins with strong evidence for longitudinal stability are less likely to fall in enhancer regions compared with DNA methylation bins with weaker evidence for longitudinal stability. The relationship with promoters, excluding the last category (where R repeated Ͼ 0.7) shows the opposite trend (Spearman's ϭ 0.94, P Ͻ .005), but the trend is weaker if the last category is included. A similar test of the lsBIN distribution in CGIs and shores reveals a significant trend concerning the islands. In CGIs the content of lsBIN increases with their longitudinal stability, from 0.17% at R repeated Ͼ 0.11 to 1.83% at R repeated Ͼ 0.70 (by Spearman's ϭ 0.945, P ϭ .00048), and this trend as expected is similar to that observed for promoter regions.

Discussion
We believe this is the first DNA methylation study to focus on SMM in a general human population. Implementing Me-DIP-seq technology, we carried out a genome-wide association analysis of DNA methylation levels with SMM variation in a dataset of 1550 middle-aged females. Due to the burden of multiple testing incorporating over 10 million tests, we restricted the methylome to the methylation signals that were minimally prone to random fluctuations. We then used these selected lsBINs in an association study with SMM in 2 independent samples, to identify robust methylation signals. press.endocrine.org/journal/jcem After QC and analysis of longitudinal stability, our data comprised of 723 029 bins showing statistically significant correlations between repeated measurements over time. The ratio of these correlations, R MZ to R DZ , also increased linearly with R repeated , suggesting that genetic factors determined this stability, as expected. This idea was further confirmed in our family-based variance component analysis of the bins most significantly associated with SMM (Supplemental Table 3), showed heritability estimates reaching 0.753 Ϯ 0.091. Thus, our data suggest that additive genetic factors affect variation of the DNA methylation levels of only small portion of methylation sites in human genome and the magnitude of this effect determines the longitudinal stability of the bins variation.
Several studies undertaken during the past decade have estimated DNA methylation twin heritability to variation in DNA methylation at individual CpG sites across the genome (eg, Refs. 29,30). Comparing the correlations in DNA methylation levels in blood samples between MZ and DZ twins' genome wide (26,690 CpG sites), Bell et al (24) estimated mean heritability of 0.18 (95% confidence interval, 0.168 -0.185). This was confirmed by Illumina HumanMethylation450 array on peripheral blood leukocytes in a family-based sample (30). The narrow sense methylation heritability in a genome-wide analysis of 4 brain regions showed a wide range of estimates (31). The mean estimate of heritability for all available loci in the latter study was only 2.8%. However, when only loci deemed heritable were included, the mean heritability estimate was 29.9%, ranging from almost 0 to almost 1.0 (31). Remarkable was the fact that the distribution pattern of their heritability estimates was virtually identical in 15 469 loci in CGIs and in 5531 loci mapped to non-CGIs. Interestingly, only about 4% of the tested loci, in this study (31) were considered heritable. The accurate estimate of genome DNA methylation depends how exactly this proportion is calculated but is of the same order of magnitude as that observed in our study.
As mentioned the higher the heritability estimate, the more longitudinally stable is methylation site (Figure 2). Furthermore, in our study the correlations between MZ twins tended to be higher for bins located in CGIs and CGI shores in comparison with "open sea" regions (Table 3), which has not been previously explored in depth due to the   limited coverage of the Illumina 450k array. In a similar vein, we observed trends with respect to the longitudinal stability of methylation bins that mapped to enhancers and promoters. The percentage of bins mapped to enhancers decreased almost monotonically with increase of R repeated , whereas promoters and CGIs displayed an opposite trend, as expected. These observations are in line with the observed epigenetic variability and stability at enhancers and promoters, but also suggest that longitudinally stable methylation in enhancers is a rare phenomenon that deserves further study. It is well established that enhancers serve as distal regulators of gene expression. Enhancer methylation is drastically altered in cancers and is closely related to altered expression profiles of cancer genes (eg, Refs. [32][33][34], suggesting that enhancer hypo-/hypermethylation could be critical in transition from normal to aberrant status of the cells. Thus, it looks that the low content of longitudinally stable methylation sites mapped to enhancers play a protective role in transcription regulation of the structural genes. Our 2-stage association analysis suggests that the methylation status of 5-7 regions is significantly associated with variation of SMM. The peak 7 regions mapped predominantly within gene bodies (4 of 7) and near to the gene transcription start site (TSS) (within 10 kb), with the exception of the intergenic signal nearest to the CAND1 gene (32 kb to TSS). Of the 7 reported genes, the topranked SMM signal was located within the DNAH12 gene, which is expressed in vastus lateralis muscle tissue after 20 weeks of endurance exercise training (35). In mice, Zfp64 has been demonstrated to participate in regulation of mesenchymal cell differentiation by suppression of myogenic, and promotion of osteoblastic, cell fate through Notch signaling (36). Tp63 has been shown to be a marker for skeletal muscle differentiation (37) and is also used for defining subsets of muscle-invasive bladder cancers (38). Finally, GUSBP11 has been listed in the C480 FANTOM5 coexpression cluster with numerous muscle-related topranked ontology terms such as skeletal muscle, contractile cell, myoblasts, striated muscle tissue, myotome, etc (39).
Hence, although the current study was conducted using whole blood, 4 of the 7 reported genes had previously been recognized in muscle-related studies, suggesting that the methylome of whole blood is informative of SMM variation. Currently, as far as we aware, no SNPs association at genome-wide association study (GWAS) significant level (P Ͻ 5 ϫ 10 Ϫ8 ) have been reported. Nevertheless, of the several recent studies, one by Guo et al (40) provides evidence of significant association of 3 SNPs (rs2507838, rs7116722, and rs11826261) in/near glycine-N-acyltransferase (GLYAT) gene in Chinese and in Caucasian populations. However, this gene is mapped, too far away (ϳ2000 kb) to consider as a reliable candidate. Thus, clearly additional GWAS and EWAS (epigenome-wide association study) are needed to understand sarcopenia.
There are certain limitations in the present study. The main potential bias concerns the usage of blood. Whole blood includes a heterogeneous group of cells whose composition can change substantially in response to a variety of intrinsic and environmental factors, eg, infection. We assessed the association between DNA methylation levels at the 7 SMM-associated and replicated DNA methylation loci with blood cell counts available for the proportion of lymphocytes, neutrophils, basophils, and eosinophils in a subset. We found no evidence that variability in WBC subtypes has a major effect on our top-ranked muscle mass related signals. The fact that we were also able to replicate these association results in a large independent sample was reassuring, and suggests that whole blood may be informative of SMM variation. MeDIP-seq provides a costeffective approach to exploring broad methylation patterns over the genome, but the resolution of the signal is not at single base pair level CpG sites. This is an additional major limitation in interpreting the methylation signals, and requires a follow up of the identified signals at greater resolution in future studies. Additional limitations include using a female population only, because males were not examined. However, using only females, could have possible advantages This approach is not prone to potential sex differences in endocrine milieu, especially sex steroids, lifestyle factors, etc, and their interactions with genotype and epigenotype, which tend to be the rule rather than the exception (41,42). Physical activity level, other lifestyle factors, medications and other environmental factors were not taken into account in this study, and may have affected the magnitude of the estimated effects. Finally, there are also some limitations related to the power of the analysis. To achieve EWAS significance accounting for multiple testing at millions of DNA methylation CpG sites would require very large sample sizes or very large effects, which have not been observed to date. To this end our association analysis was based only the longitudinally stable and most robust fraction of the methylome (Ͻ6.5% of the initial methylome).

Conclusions
This large EWAS study confirmed a major involvement of the genetic factors (heritability) to variation of SMM in a normal human population, and suggested that methylation levels of some genomic regions could underlie this variation. Four of our 7 muscle mass associated methylation signals were in or near genes previously linked to muscle, and in particular the top-ranked DNAH12 gene. Our approach based on a selection of the lsBINs reduces substantially the multiple testing problem and proved to be effective in identification of the methylation sites potentially affecting SMM variation.