Follow us on :


Take a look at the Recent articles

An epigenome-wide association study based on cell type-specific whole-genome bisulfite sequencing: Screening for DNA methylation signatures associated with bone mass

Shohei Komaki

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

E-mail : bhuvaneswari.bibleraaj@uhsm.nhs.uk

Hideki Ohmomo

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Division of Biobank and Data Management, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Tsuyoshi Hachiya

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Ryohei Furukawa

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Yuh Shiwa

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Division of Biobank and Data Management, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Mamoru Satoh

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Division of Biobank and Data Management, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Ryujin Endo

Division of Public Relations and Planning, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Division of Medical Fundamentals for Nursing, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Minoru Doita

Department of Orthopaedic Surgery, School of Medicine, Iwate Medical University, 19-1 Uchimaru, Morioka, Iwate 020-8505, Japan

Makoto Sasaki

Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Division of Ultrahigh Field MRI, Institute for Biomedical Sciences, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

Atsushi Shimizu

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University, 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

DOI: 10.15761/IMM.307

Article
Article Info
Author Info
Figures & Data

Abstract

Bone mass can change intra-individually due to aging or environmental factors. Understanding the regulation of bone metabolism by epigenetic factors, such as DNA methylation, is essential to further our understanding of bone biology and facilitate the prevention of osteoporosis. To date, a single epigenome-wide association study (EWAS) of bone density has been reported, and our knowledge of epigenetic mechanisms in bone biology is strictly limited. Here, we conducted an EWAS based on cell type-specific whole-genome bisulfite sequencing (WGBS) of CD3+/CD4+ T cells and CD14++/CD16- monocytes collected from 102 Japanese individuals, and screened DNA methylation signatures associated with bone mass. Analyses based on each cell type revealed that distinct sets of DNA methylation signatures were associated with bone mass. Some genes annotated to those DNA methylation signatures have known cell type-specific roles in bone metabolism. The results of our comprehensive screening might also contain additional novel bone-related loci, which could further our understanding of the epigenetic mechanisms of bone metabolism. With few exceptions, the cell type-specific DNA methylation signatures identified in this study are not covered by widely used arrays. Our WGBS-based EWAS highlights the importance of cell type-specific analysis with broad genome coverage, especially for discovery phase.

Key words

bone mass, CD3+/CD4+ T cell, CD14++/CD16- monocyte, epigenome-wide association study, whole genome bisulfite sequence

Introduction

The skeletal system is the most important feature supporting physical activity in human beings. Bone fracture seriously limits basic activity, and decreased bone mass is one of the major risk factors for fracture. Thus, preventing decreases in bone mass is crucial for sustaining the daily life of humans. Genetic mechanisms of osteogenesis have been well studied, and to date, various genetic markers have been reported to be associated with bone mass and density through genome-wide association studies (GWASs) [1–3]. However, bone mass and density are decreased temporarily as a result of posteriori factors, including aging, calcium and vitamin D intake, body weight, and physical activity [4–8]. Therefore, an association study focusing on epigenetic features, such as DNA methylation, which can change according to lifestyle and/or aging, should provide new insights into posteriori epigenetic mechanisms affecting bone formation and homeostasis.

DNA methylation (DNAm) is a biochemical alteration of DNA that can modulate gene transcription. Unlike nucleotide sequences, DNAm patterns vary depending on environmental factors [9]. Therefore, through the analysis of DNAm patterns, we can infer the effect of environmental factors on gene expression and subsequent phenotypes.

To date, only a single epigenetic-wide association study (EWAS) on bone mineral density has been published. In the study, the authors used a bead array (HM450, Illumina Infinium HumanMethylation450 BeadArray) and identified a single DNAm signature associated with bone density [10]. In array-based experimental systems, the number of target CpG sites is limited and increases the likelihood of overlooking potentially significant CpG sites. HM450, the most widely applied platform for human EWAS, includes approximately 480,000 probes, which cover only 1.7 % of the 28.0 million CpG sites that exist in the human genome (GRCh37) [11]. Illumina Infinium Methylation EPIC BeadChip, is a newly designed microarray-based platform, but still covers only 3.0 % of all human genome CpG sites [12]. Reduced-representation bisulfite sequencing (RRBS) covers 10% of CpG sites in the genome [13] and methyl-capture sequencing, for example SureSelect Human Methyl-Seq (Agilent Technologies), is designed to capture 4.8 million CpG sites within ± 200 kbp of target regions which cover 17.1% of all CpG sites in the genome [14]. Methylated-cytosine enrichment methods such as methylated DNA immunoprecipitation sequencing (MeDIP-Seq) and methylated DNA binding domain sequencing (MBD-Seq) can cover ~50% of all CpG sites in the human genome [13]. However, MeDIP-Seq and MBD-Seq measure the relative enrichment of methylated DNA, are less accurate for quantifying DNA methylation levels [15], and their DNAm quantifications lack single-nucleotide resolution [16]. Meanwhile, whole-genome bisulfite sequencing (WGBS) can cover ~95% of all CpG sites in the human genome with single-nucleotide resolution [17].

The bone mineral density EWAS performed by Morris et al. [10] used whole blood samples. However, they emphasized the need for a cell type-specific study because different blood cell types could have different roles in bone biology. For example, bone cell activities are regulated by cytokines derived from lymphocytes or macrophages that affect bone metabolism (osteoimmunological regulation) [18,19]. Furthermore, osteoclasts differentiate from monocytes, which means that monocytes have a direct relationship with bone cells [20,21]. Therefore, cell type-specific signatures associated with bone metabolism can be attenuated by analyzing the mixture of heterogeneous cells within whole blood samples.

Here, to increase our understanding of bone biology and the role of epigenetics in bone metabolism, we screened DNAm signatures associated with bone mass using WGBS-based EWAS in CD3+/CD4+ T cells and CD14++/CD16- monocytes isolated from 102 Japanese participants.

Materials and methods

Subjects

In this study, 109 apparently healthy subjects were enrolled. The participants are part of the Tohoku Medical Megabank Community-Based Cohort Study (TMM CommCohort Study). Demographic characteristics of the study subjects are described in Table 1. All 109 subjects were residents of the Iwate Prefecture, Japan and recruited at the Yahaba Center in the Yahaba Town, Iwate Prefecture. Details of the recruitment process were previously reported [22]. All participants provided written informed consent, and the study was approved by the Ethics Committee of Iwate Medical University (Approval ID: HG H25-19).

Bone mass measurements

A dual-energy X-ray absorptiometry (DXA) is commonly used for bone mass measurements [23]. However, the TMM CommCohort Study employed quantitative ultrasonometry (QUS) of calcaneus using an ultrasound device, Benus evo ultrasound Bone Densitometer (Nihon Kohden Co., Tokyo, Japan) to evaluate the bone mass because the equipment is of transportable size and avoids the use of ionizing radiation. Benus measures the trabecular bone area ratio that is significantly and positively associated with bone mineral density as measured by the DXA method (R2 = 0.351, P < 0.001) [24]. Therefore, in this study, the bone area ratio was used to represent bone mass. For 109 subjects, the trabecular bone mass was measured by an identical ultrasound device at the Yahaba Center following the manufacturer’s instructions.

Sample preparation and DNA methylation profiling

Detailed methods for sample preparation and DNAm profiling were previously described by Hachiya et al. [25]. In brief, for the isolation of peripheral blood mononuclear cells, 8 ml of whole blood was collected from each participant, and CD3+/CD4+ T cells (CD4T) and CD14++/CD16- monocytes were further isolated. Genomic DNA extracted from each cell type went through bisulfite conversion and library preparation, followed by WGBS. After removal of the adapter, sequence reads (length ≥ 20 bp) were mapped onto the GRCh37d5 reference genome. Overlaps between paired-end reads were clipped. Then, DNAm levels at each CpG site were estimated by dividing the number of the mapped reads that contained unconverted cytosine by the total number of mapped reads. CpG sites with extreme depth (< 6 or > 300) or with a high frequency of missing data among subjects (call rate < 50%) were excluded from the dataset.

Epigenome wide association analysis

We performed six epigenome-wide association analyses based on DNAm profiles of: (1) CD4T of all subjects, (2) CD4T of females, (3) CD4T of males, (4) monocytes of all subjects, (5) monocytes of females, and (6) monocytes of males. A linear regression (LM) model was used to test the relationship between the DNAm level of each ~24 million CpG sites and bone mass. CpG sites on sex chromosomes and mitochondrion were excluded from the analyses. The LM was performed using the lm function of the R version 3.2.0, specifying bone mass as the dependent variable and DNAm level, age, and BMI as the explanatory variables. Analyses (1) and (4) were further adjusted by sex. Each analysis was assessed based on the genomic inflation factor (λ) and quantile-quantile plot (QQ plot). In this study, a less stringent threshold of P value = 1.0 × 10-7, suggested for the discovery criteria, was applied [26].

Results

Mean and standard deviation (SD) of bone mass (%) among all 109 subjects were 28.94 and 4.51, respectively (Table 1). Females showed lower bone mass mean value and larger sd than did males (Table 1) (F test P value = 0.032, Wilcoxon rank sum test P value = 0.009). The bone mass of females greatly fluctuates in association with estrogen level, especially at the early menopause stage [27,28], which can explain the observed gender differences in bone mass. These results imply that there are sex-specific mechanisms characterizing the bone mass of study subjects, and emphasize the need for sex-specific analyses to explore the epigenetic regulation of bone metabolism.

Table 1. Description of the cohort studied

 

CD4T

Monocytes

N (% Female)

102 (52.90)

102 (53.90)

Bone area ratio (%) (mean ± SD)

28.94 ± 4.59

28.99 ± 4.58

Age (year) (mean ± SD)

59.01 ± 11.30

59.50 ± 10.87

Number of CpG sites analyzed

24,036,660

23,940,752

From 109 subjects, 102 and 102 were used for the analyses of CD4T and monocytes, respectively, with 95 subjects overlapping. SD: standard deviation, CD4T: CD3+/CD4+ T cells.

From the 109 subjects, DNAm profiles of CD4T and monocytes were obtained from 102 and 102 subjects, respectively (Table 1), with 95 subjects overlapped. For CD4T, 24,036,660 CpG sites were analyzed, and 23,940,752 sites were analyzed for monocytes. Resultant P values of each analysis, based on different datasets, are presented as Manhattan and QQ plots (Figure 1 and Supplementary Figure S1). Genomic inflation factors for each analysis approximate to 1, indicating that confounding factors were well controlled in the analyses (Supplementary Figure S1).

Figure 1. Manhattan plots for each association analysis based on different cell types and datasets. The red line indicates P value = 1.0 × 10-7. CD4T: CD3+/CD4+ T cells.

CpG sites with P values below the threshold of 1.0 × 10-7 in each analysis are listed in Tables 2 and 3. In analyses based on CD4T of both sexes, females only, and males only, 18, 4, and 9 CpG sites showed P values < 1.0 × 10-7, respectively. Analyses based on monocytes of both sexes, females only, and males only identified 18, 22, and 7 CpG sites with P values < 1.0 × 10-7, respectively.

Table 2. CpG sites associated with bone mass (P < 1.0 × 10-7) identified in CD3+/CD4+ T cells

Dataset

Position*1

DNAm level (mean ± SD)

Coefficient (95% CI)

P value

Gene

Both sexes

12:100925094

75.80 ± 29.19

-3.90 (-5.07 – -2.73)

2.15 × 10-9

NR1H4

 

16:69898597

97.71 ± 5.11

-0.71 (-0.93 – -0.50)

2.85 × 10-9

WWP2*2

 

20:2119971

97.59 ± 2.92

-0.38 (-0.50 – -0.26)

1.22 × 10-8

STK35

 

22:49949200

95.19 ± 5.08

-0.64 (-0.85 – -0.43)

1.69 × 10-8

RP1-29C18.10/C22orf34

 

9:34115355

94.55 ± 6.92

-0.89 (-1.17 – -0.60)

1.73 × 10-8

DCAF12

 

1:78131263

98.89 ± 3.26

-0.46 (-0.61 – -0.32)

1.79 × 10-8

ZZZ3

 

3:54141925

97.46 ± 7.70

-1.14 (-1.50 – -0.79)

2.21 × 10-8

 

 

11:10061108

97.46 ± 4.62

-0.56 (-0.74 – -0.38)

2.61 × 10-8

SBF2*3

 

1:246227233

97.04 ± 3.09

-0.39 (-0.52 – -0.27)

2.62 × 10-8

SMYD3

 

17:980028

98.17 ± 2.88

-0.36 (-0.48 – -0.24)

3.09 × 10-8

ABR

 

13:98739376

92.83 ± 9.98

-1.22 (-1.63 – -0.81)

4.52 × 10-8

FARP1

 

22:49831298

92.28 ± 4.83

-0.60 (-0.80 – -0.40)

4.74 × 10-8

C22orf34

 

22:49967265

93.68 ± 6.72

-0.83 (-1.11 – -0.55)

4.90 × 10-8

RP1-29C18.9/C22orf34

 

1:145549040

1.07 ± 2.81

0.35 (0.23 – 0.46)

4.99 × 10-8

ANKRD35

 

11:609352

96.70 ± 6.47

-0.80 (-1.07 – -0.53)

6.28 × 10-8

PHRF1

 

5:176738773

0.47 ± 1.07

0.13 (0.09 – 0.18)

7.26 × 10-8

MXD3

 

3:121842023

96.57 ± 3.27

-0.40 (-0.54 – -0.27)

7.56 × 10-8

CD86*3

 

2:226290206

95.65 ± 6.08

-0.75 (-1.01 – -0.49)

9.82 × 10-8

NYAP2

Female

1:21478594

95.61 ± 7.46

-1.38 (-1.78 – -0.98)

1.18 × 10-8

EIF4G3

 

19:10127122

98.74 ± 2.59

-0.46 (-0.60 – -0.32)

1.80 × 10-8

RDH8

 

5:77319949

98.40 ± 4.02

-0.68 (-0.89 – -0.48)

2.15 × 10-8

AP3B1

 

4:68576716

99.05 ± 2.77

-0.49 (-0.64 – -0.33)

6.60 × 10-8

UBA6-AS1

Male

3:197772225

96.91 ± 3.43

-0.70 (-0.86 – -0.53)

9.36 × 10-11

LMLN

 

5:99290935

96.48 ± 2.76

-0.50 (-0.64 – -0.35)

9.22 × 10-9

CTD-2160D9.1

 

5:26240956

97.78 ± 3.46

-0.63 (-0.82 – -0.44)

2.59 × 10-8

RP11-351N6.1

 

2:184159598

95.89 ± 3.58

-0.65 (-0.84 – -0.45)

3.80 × 10-8

LIN28AP1

 

19:8455454

0.61 ± 1.52

0.27 (0.19 – 0.36)

4.92 × 10-8

RAB11B/RAB11B-AS1

 

9:126673994

96.62 ± 5.29

-0.99 (-1.29 – -0.68)

5.63 × 10-8

DENND1A

 

8:38854608

0.45 ± 1.40

0.25 (0.17 – 0.33)

6.76 × 10-8

ADAM9*3

 

12:69611102

94.12 ± 5.56

-0.98 (-1.29 – -0.67)

8.77 × 10-8

RP11-324P9.1

 

2:110449774

84.72 ± 9.25

-1.66 (-2.19 – -1.13)

9.34 × 10-8

BMS1P19

*1 Chromosome number and position of the CpG site.

*2 Gene with evident roles in bone development and/or homeostasis.

*3 Gene with a suggested function in, or relationship with, bone development and/or homeostasis.

DNAm: DNA methylation, SD: standard deviation, CI: confidence interval.

Table 3. CpG sites associated with bone mass (P < 1.0 × 10-7) identified in monocytes

Dataset

Position*1

DNAm level (mean ± SD)

Coefficient

(95% CI)

P value

Gene

Both sexes

2:61108446

0.29 ± 1.38

0.21 (0.15 – 0.26)

2.91 × 10-11

AC010733.4

 

22:47133916

0.84 ± 2.80

0.40 (0.28 – 0.51)

3.73 × 10-10

CERK

 

2:121815071

96.54 ± 4.01

-0.56 (-0.72 – -0.39)

5.95 × 10-10

Y RNA

 

3:188489102

98.19 ± 3.52

-0.46 (-0.60 – -0.32)

3.02 × 10-9

LPP

 

8:8750676

1.11 ± 1.71

0.22 (0.15 – 0.30)

7.51 × 10-9

MFHAS1

 

14:36295936

0.51 ± 1.43

0.19 (0.13 – 0.25)

1.30 × 10-8

RP11-317N8.5/BRMS1L

 

2:195259677

85.29 ± 14.8

-2.22 (-2.89 – -1.55)

2.02 × 10-8

AC018799.1

 

2:69941555

97.83 ± 3.04

-0.39 (-0.52 – -0.26)

2.66 × 10-8

ANXA4

 

3:173900045

85.08 ± 18.96

-2.76 (-3.64 – -1.89)

2.95 × 10-8

NLGN1

 

2:190657692

96.39 ± 3.84

-0.48 (-0.64 – -0.32)

4.33 × 10-8

PMS1

 

12:56474020

0.25 ± 1.11

0.14 (0.09 – 0.18)

6.30 × 10-8

ERBB3*3

 

17:1535968

92.9 ± 13.44

-1.63 (-2.19 – -1.08)

6.43 × 10-8

SCARF1

 

4:99080231

92.82 ± 9.63

-1.41 (-1.86 – -0.95)

6.45 × 10-8

Y RNA

 

8:21530651

94.44 ± 12.01

-1.60 (-2.13 – -1.07)

6.62 × 10-8

GFRA2

 

2:240322590

0.50 ± 1.79

0.22 (0.15 – 0.30)

6.72 × 10-8

HDAC4*3

 

12:129772248

90.49 ± 7.00

-0.85 (-1.14 – -0.56)

8.41 × 10-8

TMEM132D

 

3:16191417

96.63 ± 7.87

-1.00 (-1.34 – -0.66)

9.09 × 10-8

GALNT15

 

14:39639583

1.60 ± 2.61

0.32 (0.21 – 0.43)

9.11 × 10-8

TRAPPC6B

Female

6:35995637

0.37 ± 1.07

0.21 (0.16 – 0.25)

2.44× 10-11

MAPK14*2

 

18:59560980

0.51 ± 1.48

0.27 (0.20 – 0.35)

1.87 × 10-9

RNF152

 

6:163371471

94.89 ± 8.66

-2.42 (-2.96 – -1.87)

2.54 × 10-9

PACRG

 

21:17220599

97.23 ± 4.77

-0.89 (-1.14 – -0.63)

5.16 × 10-9

USP25

 

10:60936039

6.68 ± 5.35

0.98 (0.70 – 1.26)

6.22 × 10-9

PHYHIPL

 

11:497792

96.88 ± 4.27

-0.78 (-1.01 – -0.56)

9.53 × 10-9

RNH1*3

 

4:81105396

0.52 ± 1.96

0.35 (0.25 – 0.45)

1.84 × 10-8

RP11-377G16.2/PRDM8

 

5:65423122

84.35 ± 21.30

-3.75 (-4.89 – -2.60)

3.07 × 10-8

SREK1

 

4:76766071

97.61 ± 3.13

-0.55 (-0.72 – -0.38)

3.34 × 10-8

RP11-556N4.1

 

16:56225078

0.77 ± 2.84

0.50 (0.36 – 0.65)

3.68 × 10-8

RP11-461O7.1

 

20:13765582

0.34 ± 1.48

0.26 (0.18 – 0.34)

3.78 × 10-8

NDUFAF5

 

1:201449877

0.23 ± 0.96

0.17 (0.12 – 0.22)

3.86 × 10-8

CSRP1

 

11:124035292

98.58 ± 3.03

-0.53 (-0.70 – -0.37)

4.21 × 10-8

OR10D1P

 

16:15188395 (cg16603896)

0.96 ± 3.22

0.54 (0.37 – 0.71)

4.60 × 10-8

RP11-72I8.1/PDXDC1

 

2:24660027

92.65 ± 8.78

-1.53 (-2.00 – -1.06)

4.83 × 10-8

NCOA1

 

5:64919917

0.23 ± 0.96

0.16 (0.11 – 0.21)

6.02 × 10-8

TRIM23

 

6:36954391

0.47 ± 1.62

0.29 (0.20 – 0.39)

6.79 × 10-8

MTCH1

 

7:29157080

98.58 ± 3.53

-0.61 (-0.81 – -0.42)

7.31 × 10-8

CPVL

 

2:10338969

97.26 ± 3.76

-0.64 (-0.85 – -0.44)

8.54 × 10-8

C2orf48

 

4:4377945

95.56 ± 3.49

-0.61 (-0.80 – -0.41)

8.58 × 10-8

NSG1

 

1:23857808

0.97 ± 3.41

0.72 (0.50 – 0.93)

9.08 × 10-8

E2F2

 

2:230125683

84.64 ± 9.88

-1.72 (-2.27 – -1.17)

9.61 × 10-8

PID1

Male

10:60540166

87.39 ± 10.86

-2.06 (-2.64 – -1.48)

6.42 × 10-9

BICC1*2

 

1:3463645

93.37 ± 4.07

-0.75 (-0.96 – -0.54)

6.84 × 10-9

MEGF6

 

19:47039387

98.04 ± 4.18

-0.84 (-1.08 – -0.59)

3.90 × 10-8

PPP5D1

 

2:172721808

91.63 ± 11.99

-2.36 (-3.07 – -1.65)

4.68 × 10-8

SLC25A12

 

12:6588863

96.45 ± 3.77

-0.66 (-0.86 – -0.46)

5.54 × 10-8

RP1-102E24.6

 

20:50007014

96.99 ± 4.16

-0.73 (-0.96 – -0.50)

7.50 × 10-8

NFATC2*3

 

14:93054727

92.29 ± 19.21

-3.50 (-4.60 – -2.40)

8.00 × 10-8

RIN3*3

*1 Chromosome number and the position of CpG site.

*2 Gene with evident roles in bone development and/or homeostasis.

*3 Gene with suggested function in, or relationship with, bone development and/or homeostasis. DNAm: DNA methylation, SD: standard deviat

Discussion

We screened CpG signatures associated with bone mass by means of WGBS-based EWAS. As a result, different sets of CpG signatures were identified among cell type- and sex-specific analyses. Genes located close (within ± 1 Mbp) to CpG sites with P values < 1.0 × 10-7 are listed in Tables 2 and 3. Among them, WWP2, BICC1, and MAPK14, have known biological functions related to bone development or homeostasis. WWP2 encodes a member of the NEDD4 family of E3 ligases that function in protein ubiquitylation. Indeed, monoubiquitylation of goosecoid by WWP2 protein in chondrocytes regulates craniofacial skeletal patterning [29]. Furthermore, knockdown and overexpression experiments of Wwp2 in mesenchymal stem cells of mice by Zhu et al. [30] suggested that WWP2 behaves as positive regulator for osteogenesis through ubiquitylation of RUNX2, which is required for the early stage of osteoblast differentiation. Therefore, suppressed transcription of WWP2 can be responsible for lowered bone mass, and might explain the negative association observed between bone mass and WWP2 DNAm level in CD4T.

BICC1 is a gene encoding an RNA-binding protein. Knockdown experiments in mice and GWAS of human bone mineral density show that Bicc1/BICC1 plays a significant role in osteoblastogenesis and affects bone mineral density [31]. Briefly, BICC1 positively regulates osteoblastogenesis, and suppression of BICC1 expression leads to decreased bone mineral density. In both mice and humans, the association between Bicc1/BICC1 and bone mineral density was demonstrated in male subjects. Consistently, a significant association was observed between BICC1 DNAm and bone mass in monocytes of male subjects in this study. The negative association between DNAm level of the CpG on BICC1 and bone mass implies that methylation of the site decreases BICC1 expression, leading to lowered bone mass in male subjects. Of all blood cell types, only plasma and monocyte show expression of BICC1 (GeneCards database, http://www.genecards.org; accessed on July 1, 2017), which could explain why a BICC1 DNAm signature was identified from monocytes but not from CD4T cells in this study.

Mitogen-activated protein kinase 14 (p38α), encoded by MAPK14, has a function in bone resorbing. Expression of genes essential for osteoclast differentiation, such as SP7, ALPL, and BGLAP, is regulated by p38α [32]. Using a specific p38α kinase inhibitor, osteoclast differentiation is suppressed. Furthermore, it is suggested that phosphorylation of microphthalmia-associated transcription factor, required for the osteoclast maturation, is stimulated by phosphorylated p38α [33]. Since osteoclasts function in bone resorption, stimulated osteoclast formation can lead to excess bone loss. The function of p38α to negatively affect osteogenesis through activating osteoclastogenesis and osteoclast function might explain the negative association between bone mass and CpG DNAm levels of MAPK14. This association was observed only from the analysis based on monocytes of females. This is consistent with the facts that osteoclasts are derived from monocytes [20,21], and that the p38α pathway is an important regulator of osteoporosis in postmenopausal women and induced by estrogen deficiency [28].

Other genes listed in Table 2 are thought to function in bone metabolism, or are associated with bone density through GWAS (ERBB3, [34]; SBF2, [35]; RNH1, [36]; NFATC2, [37]; RIN3, [38]; ADAM9, [39]; HDAC4, [40]; and CD86, [41]).

Furthermore, genes other than above mentioned may have unreported functions in bone metabolism including osteogenesis and osteoclastogenesis. For example, NR1H4, also known as farnesoid X receptor alpha (FXRα), is a member of the nuclear receptor superfamily of ligand-dependent transcription factors. A CpG site, chr12:100925094, whose DNAm level is significantly associated with bone mass, is located on NR1H4. This receptor is an activator of gluconeogenic pathways and a known regulator of glucocorticoid receptor expression and activity [42]. Glucocorticoids have a detrimental effect on bone formation and homeostasis through suppressing the differentiation and replication of osteoblast cells, and by inhibiting the function of mature osteoblast cells [43]. In addition, increased osteoblast apoptosis is induced by glucocorticoids and glucocorticoids decrease the function of osteocytes [43]. Given their anti-inflammatory and immunosuppressive properties, glucocorticoids have been widely used as therapeutic agents, even though glucocorticoid therapy can cause bone loss and fractures in patients. The negative correlation between bone mass and DNAm level of NR1H4 CpG observed in this study implies that DNAm suppresses the transcription or function of NR1H4, resulting in bone loss.

A previous EWAS of bone mineral density reported eight CpG sites as suggestive DNAm signatures in individual cohorts but did not identify any DNAm signatures common across the studied cohorts [10]. In this study, these eight previously reported CpG sites showed no significant associations with bone mass (P values of > 0.130) (Supplementary Table S1). However, we analyzed Japanese individuals applying bone area ratio of calcaneus measured by QUS as bone mass, and DNAm level was measured using purified CD4T and monocyte cells. The different materials and methods used in the two studies might have been the cause for the differing results. Thus, our results do not deny possible relationships between bone mineral density and DNAm levels of the CpG sites suggested by Morris et al. [10].

We note that the number of subjects enrolled in this study was around 100, and trabecular bone area ratio of calcaneus measured by QUS was applied as a phenotypic variable. Therefore, our study was based on limited data. Furthermore, the association analyses were performed based on a single Japanese cohort and associations between DNAm levels and bone mass were not validated in any other independent cohorts. Therefore, associations reported in this study remain to be confirmed by additional studies.

Despite the above-mentioned limitations, whole-genome DNAm analyses applied in this study might provide new insights into the role of DNAm in bone metabolism. Among 78 DNAm signatures identified in the present EWAS, only a single CpG site [chr16:15188395 (cg16603896); Table 3] is a target of the most commonly applied platform for DNA methylation studies in humans, HM450. In this study, furthermore, several CpG sites were identified as DNAm signatures with small P values, while neighboring CpG sites surrounding the signatures did not show such associations between DNAm level and bone mass (Supplementary Figure S2). This result implies that, without the high genome-coverage, DNAm association studies could have overlooked a large number of potential DNAm signals. Therefore, this study emphasizes the importance of WGBS-based analysis with the highest genome-coverage, especially for discovery phrase requiring comprehensive screening with few genomic region omissions.

The aim of this study was to screen DNA methylation signatures associated with bone mass using cell type-specific analyses with the widest genome coverage. We have successfully conducted analyses and have presented candidate DNAm signatures. EWAS with larger-scale cohorts focusing on the specific regions identified in this study will further elucidate the role of epigenetic mechanisms in bone biology.

Acknowledgements

The authors appreciate the participants of the Tohoku Medical Megabank Project. We would like to thank the members of the Iwate Tohoku Medical Megabank Organization of Iwate Medical University (IMM) and the Tohoku Medical Megabank Organization of Tohoku University (ToMMo) for their encouragement and support. We deeply thank Ms. Kanako Ono for her assistance with sample processing and handling.

Funding information

  • The Reconstruction Agency, the Ministry of Education, Culture, Sports, Science and Technology (MEXT)
  • The Japan Agency for Medical Research and Development (AMED)

Competing interest

None

View supplementary data

References

  1. Cho YS, Go MJ, Kim YJ, Heo JY, Oh JH, et al. (2009) A large-scale genome-wide association study of Asian populations uncovers genetic factors influencing eight quantitative traits. Nat Genet 41: 527–534. [Crossref]
  2. Rivadeneira F, Styrkársdottir U, Estrada K, Halldórsson BV, Hsu YH, et al. (2009) Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. Nat Genet 41: 1199–1206. [Crossref]
  3. Estrada K, Styrkarsdottir U, Evangelou E, Hsu YH, Duncan EL, et al. (2012) Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture. Nat Genet 44: 491-501. [Crossref] 
  4. Demontiero O, Vidal C, Duque G (2012) Aging and bone loss: new insights for the clinician. Ther Adv Musculoskelet Dis 4: 61-76. [Crossref] 
  5. Dawson-Hughes B, Harris SS, Krall EA, Dallal GE (1997) Effect of calcium and vitamin D supplementation on bone density in men and women 65 years of age or older. N Engl J Med 337: 670-676. [Crossref] 
  6. Booth SL, Broe KE, Gagnon DR, Tucker KL, Hannan MT, et al. (2003) Vitamin K intake and bone mineral density in women and men. Am J Clin Nutr 77: 512-516. [Crossref] 
  7. Shapses SA, Riedt CS (2006) Bone, body weight, and weight reduction: what are the concerns? J Nutr 136: 1453-1456. [Crossref] 
  8. Warden SJ, Mantila Roosa SM, Kersh ME, Hurd AL, Fleisig GS, et al. (2014) Physical activity when young provides lifelong benefits to cortical bone size and strength in men. Proc Natl Acad Sci U S A 111: 5337–5342.
  9. Foley DL, Craig JM, Morley R, Olsson CA, Dwyer T, et al. (2009) Prospects for epigenetic epidemiology. Am J Epidemiol 169: 389-400. [Crossref] 
  10. Morris JA, Tsai PC, Joehanes R, Zheng J, Trajanoska K, et al. (2017) Epigenome-wide Association of DNA Methylation in Whole Blood With Bone Mineral Density. J Bone Miner Res 32: 1644-1650. [Crossref]
  11. Bibikova M, Barnes B, Tsan C, Ho V, Klotzle B, et al. (2011) High density DNA methylation array with single CpG site resolution. Genomics 98: 288-295. [Crossref] 
  12. Pidsley R, Zotenko E, Peters TJ, Lawrence MG, Risbridger GP, et al. (2016) Critical evaluation of the Illumina MethylationEPIC BeadChip microarray for whole-genome DNA methylation profiling. Genome Biol 17: 208. [Crossref]
  13. Walker DL, Bhagwate AV, Baheti S, Smalley RL, Hilker CA, et al. (2015) DNA methylation profiling: comparison of genome-wide sequencing methods and the Infinium Human Methylation 450 Bead Chip. Epigenomics 7: 1–16. [Crossref]
  14. Teh AL, Pan H, Lin X, Lim YI, Patro CPK, et al. (2016) Comparison of Methyl-capture Sequencing vs. Infinium 450K methylation array for methylome analysis in clinical samples. Epigenetics 11: 36–48. [Crossref]
  15. Bock C, Halbritter F, Carmona FJ, Tierling S, Datlinger P, et al. (2016) Quantitative comparison of DNA methylation assays for biomarker development and clinical applications. Nat Biotechnol 34: 726–737. [Crossref]
  16. Harris RA, Wang T, Coarfa C, Nagarajan RP, Hong C, et al. (2010) Comparison of sequencing-based methods to profile DNA methylation and identification of monoallelic epigenetic modifications. Nat Biotechnol 28: 1097–1105.[Crossref]
  17. Stirzaker C, Taberlay PC, Statham AL, Clark SJ (2014) Mining cancer methylomes: prospects and challenges. Trends Genet 30: 75-84. [Crossref] 
  18. Walsh MC, Kim N, Kadono Y, Rho J, Lee SY, et al. (2006) Osteoimmunology: Interplay between the immune system and bone metabolism. Annu Rev Immunol 24: 33–63. [Crossref]
  19. Mori G, D'Amelio P, Faccio R, Brunetti G (2013) The Interplay between the bone and the immune system. Clin Dev Immunol 2013: 720504. [Crossref] 
  20. Udagawa N, Takahashi N, Akatsu T, Tanaka H, Sasaki T, et al. (1990) Origin of osteoclasts: mature monocytes and macrophages are capable of differentiating into osteoclasts under a suitable microenvironment prepared by bone marrow-derived stromal cells. Proc Natl Acad Sci U S A 87: 7260–7264. [Crossref]
  21. Boyle WJ, Simonet WS, Lacey DL (2003) Osteoclast differentiation and activation. Nature 423: 337-342. [Crossref] 
  22. Kuriyama S, Yaegashi N, Nagami F, Arai T, Kawaguchi Y, et al. (2016) The Tohoku Medical Megabank Project: Design and Mission. J Epidemiol 26: 493-511. [Crossref] 
  23. Pisani P, Renna MD, Conversano F, Casciaro E, Muratore M, et al. (2013) Screening and early diagnosis of osteoporosis through X-ray and ultrasound based techniques. World J Radiol 5: 398–410.
  24. Tatebe T, Takanashi Y, Mizutani Y, Miura T (1998) A study for validity of ultrasound bone densitometry using “Benus.” Res J Heal Sport Sci Chukyo Univ 40: 103–110.
  25. Hachiya T, Furukawa R, Shiwa Y, Ohmomo H, Ono K, et al. (2017) Genome-wide identification of inter-individually variable DNA methylation sites improves the efficacy of epigenetic association studies. npj Genomic Med 2: 11.
  26. Panagiotou OA, Ioannidis JPA, Hirschhorn JN, Abecasis GR, Frayling TM, et al. (2012) What should the genome-wide significance threshold be? Empirical replication of borderline genetic associations. Int J Epidemiol 41: 273–286.
  27. Clarke BL, Khosla S (2009) Androgens and bone. Steroids 74: 296-305. [Crossref] 
  28. Caverzasio J, Higgins L, Ammann P (2008) Prevention of trabecular bone loss induced by estrogen deficiency by a selective p38alpha inhibitor. J Bone Miner Res 23: 1389-1397. [Crossref] 
  29. Zou W, Chen X, Shim J, Huang Z, Brady N, et al. (2011) The E3 ubiquitin ligase Wwp2 regulates craniofacial development through mono-ubiquitylation of Goosecoid. Nat Cell Biol 13: 59–65. [Crossref]
  30. Zhu W, He X, Hua Y, Li Q, Wang J, et al. (2017) The E3 Uubiquitin ligase WWP2 facilitates RUNX2 transactivation through a mono-ubiquitination manner during osteogenic differentiation. J Biol Chem 292: 11178–11188. [Crossref]
  31. Mesner LD, Ray B, Hsu YH, Manichaikul A, Lum E, et al. (2014) Bicc1 is a genetic determinant of osteoblastogenesis and bone mineral density. J Clin Invest 124: 2736-2749. [Crossref] 
  32. Thouverey C, Caverzasio J (2015) Focus on the p38 MAPK signaling pathway in bone development and maintenance. Bonekey Rep 4: 711. [Crossref] 
  33. Mansky KC, Sankar U, Han J, Ostrowski MC (2002) Microphthalmia transcription factor is a target of the p38 MAPK pathway in response to receptor activator of NF-kappa B ligand signaling. J Biol Chem 277: 11077–10083. [Crossref]
  34. Lin SH, Cheng CJ, Lee YC, Ye X, Tsai WW, et al. (2008) A 45-kDa ErbB3 secreted by prostate cancer cells promotes bone formation. Oncogene 27: 5195-5203. [Crossref] 
  35. Lei SF, Tan LJ, Liu XG, Wang L, Yan H, et al. (2009) Genome-wide association study identifies two novel loci containing FLNB and SBF2 genes underlying stature variation. Hum Mol Genet 18: 1661–1669. [Crossref]
  36. Kantawong F, Burchmore R, Gadegaard N, Oreffo ROC, Dalby MJ (2009) Proteomic analysis of human osteoprogenitor response to disordered nanotopography. J R Soc Interface 6: 1075–1086. [Crossref]
  37. Sitara D, Aliprantis AO (2010) Transcriptional regulation of bone and joint remodeling by NFAT. Immunol Rev 233: 286-300. [Crossref] 
  38. Kemp JP, Medina-Gomez C, Estrada K, St Pourcain B, Heppe DHM, et al. (2014) Phenotypic dissection of bone mineral density reveals skeletal site specificity and facilitates the identification of novel loci in the genetic regulation of bone mass attainment. PLoS Genet 10: e1004423. [Crossref]
  39. Karadag A, Zhou M, Croucher PI (2006) ADAM-9 (MDC-9/meltrin-γ), a member of the a disintegrin and metalloproteinase family, regulates myeloma-cell-induced interleukin-6 production in osteoblasts by direct interaction with the avß5 integrin. Blood 107: 3271–3278. [Crossref]
  40. Vega RB, Matsuda K, Oh J, Barbosa AC, Yang X, et al. (2004) Histone deacetylase 4 controls chondrocyte hypertrophy during skeletogenesis. Cell 119: 555–566. [Crossref]
  41. Bozec A, Zaiss MM, Kagwiria R, Voll R, Rauh M, et al. (2014) T cell costimulation molecules CD80/86 inhibit osteoclast differentiation by inducing the IDO/tryptophan pathway. Sci Transl Med 6: 235–246. [Crossref]
  42. Renga B, Mencarelli A, D’Amore C, Cipriani S, Baldelli F, et al. (2012) Glucocorticoid receptor mediates the gluconeogenic activity of the farnesoid X receptor in the fasting condition. FASEB J 26: 3021–3031. [Crossref]
  43. Canalis E, Mazziotti G, Giustina A, Bilezikian JP (2007) Glucocorticoid-induced osteoporosis: pathophysiology and therapy. Osteoporos Int 18: 1319–28. [Crossref]
  44. Weisenberger DJ (2014) Characterizing DNA methylation alterations from The Cancer Genome Atlas. J Clin Invest 124: 17-23. [Crossref] 

Editorial Information

Editor-in-Chief

Ivan Gout
University College London

Ricardo H. Alvarez
Cancer Treatment Centers of America

Article Type

Research Article

Publication history

Received date: September 20, 2017
Accepted date: October 09, 2017
Published date: October 12, 2017

Copyright

© 2017 Komaki S. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Citation

Komaki S, Ohmomo H, Hachiya T, Furukawa R, Shiwa Y, et al. (2017) An epigenome-wide association study based on cell type-specific whole-genome bisulfite sequencing: Screening for DNA methylation signatures associated with bone mass. Integr Mol Med 4: DOI: 10.15761/IMM.307

Corresponding author

Atsushi Shimizu, Ph.D

Division of Biomedical Information Analysis, Iwate Tohoku Medical Megabank Organization, Disaster Reconstruction Center, Iwate Medical University 2-1-1 Nishitokuta, Yahaba, Shiwa, Iwate 028-3694, Japan

E-mail : bhuvaneswari.bibleraaj@uhsm.nhs.uk

Figure 1. Manhattan plots for each association analysis based on different cell types and datasets. The red line indicates P value = 1.0 × 10-7. CD4T: CD3+/CD4+ T cells.

Table 1. Description of the cohort studied

 

CD4T

Monocytes

N (% Female)

102 (52.90)

102 (53.90)

Bone area ratio (%) (mean ± SD)

28.94 ± 4.59

28.99 ± 4.58

Age (year) (mean ± SD)

59.01 ± 11.30

59.50 ± 10.87

Number of CpG sites analyzed

24,036,660

23,940,752

From 109 subjects, 102 and 102 were used for the analyses of CD4T and monocytes, respectively, with 95 subjects overlapping. SD: standard deviation, CD4T: CD3+/CD4+ T cells.

Table 2. CpG sites associated with bone mass (P < 1.0 × 10-7) identified in CD3+/CD4+ T cells

Dataset

Position*1

DNAm level (mean ± SD)

Coefficient (95% CI)

P value

Gene

Both sexes

12:100925094

75.80 ± 29.19

-3.90 (-5.07 – -2.73)

2.15 × 10-9

NR1H4

 

16:69898597

97.71 ± 5.11

-0.71 (-0.93 – -0.50)

2.85 × 10-9

WWP2*2

 

20:2119971

97.59 ± 2.92

-0.38 (-0.50 – -0.26)

1.22 × 10-8

STK35

 

22:49949200

95.19 ± 5.08

-0.64 (-0.85 – -0.43)

1.69 × 10-8

RP1-29C18.10/C22orf34

 

9:34115355

94.55 ± 6.92

-0.89 (-1.17 – -0.60)

1.73 × 10-8

DCAF12

 

1:78131263

98.89 ± 3.26

-0.46 (-0.61 – -0.32)

1.79 × 10-8

ZZZ3

 

3:54141925

97.46 ± 7.70

-1.14 (-1.50 – -0.79)

2.21 × 10-8

 

 

11:10061108

97.46 ± 4.62

-0.56 (-0.74 – -0.38)

2.61 × 10-8

SBF2*3

 

1:246227233

97.04 ± 3.09

-0.39 (-0.52 – -0.27)

2.62 × 10-8

SMYD3

 

17:980028

98.17 ± 2.88

-0.36 (-0.48 – -0.24)

3.09 × 10-8

ABR

 

13:98739376

92.83 ± 9.98

-1.22 (-1.63 – -0.81)

4.52 × 10-8

FARP1

 

22:49831298

92.28 ± 4.83

-0.60 (-0.80 – -0.40)

4.74 × 10-8

C22orf34

 

22:49967265

93.68 ± 6.72

-0.83 (-1.11 – -0.55)

4.90 × 10-8

RP1-29C18.9/C22orf34

 

1:145549040

1.07 ± 2.81

0.35 (0.23 – 0.46)

4.99 × 10-8

ANKRD35

 

11:609352

96.70 ± 6.47

-0.80 (-1.07 – -0.53)

6.28 × 10-8

PHRF1

 

5:176738773

0.47 ± 1.07

0.13 (0.09 – 0.18)

7.26 × 10-8

MXD3

 

3:121842023

96.57 ± 3.27

-0.40 (-0.54 – -0.27)

7.56 × 10-8

CD86*3

 

2:226290206

95.65 ± 6.08

-0.75 (-1.01 – -0.49)

9.82 × 10-8

NYAP2

Female

1:21478594

95.61 ± 7.46

-1.38 (-1.78 – -0.98)

1.18 × 10-8

EIF4G3

 

19:10127122

98.74 ± 2.59

-0.46 (-0.60 – -0.32)

1.80 × 10-8

RDH8

 

5:77319949

98.40 ± 4.02

-0.68 (-0.89 – -0.48)

2.15 × 10-8

AP3B1

 

4:68576716

99.05 ± 2.77

-0.49 (-0.64 – -0.33)

6.60 × 10-8

UBA6-AS1

Male

3:197772225

96.91 ± 3.43

-0.70 (-0.86 – -0.53)

9.36 × 10-11

LMLN

 

5:99290935

96.48 ± 2.76

-0.50 (-0.64 – -0.35)

9.22 × 10-9

CTD-2160D9.1

 

5:26240956

97.78 ± 3.46

-0.63 (-0.82 – -0.44)

2.59 × 10-8

RP11-351N6.1

 

2:184159598

95.89 ± 3.58

-0.65 (-0.84 – -0.45)

3.80 × 10-8

LIN28AP1

 

19:8455454

0.61 ± 1.52

0.27 (0.19 – 0.36)

4.92 × 10-8

RAB11B/RAB11B-AS1

 

9:126673994

96.62 ± 5.29

-0.99 (-1.29 – -0.68)

5.63 × 10-8

DENND1A

 

8:38854608

0.45 ± 1.40

0.25 (0.17 – 0.33)

6.76 × 10-8

ADAM9*3

 

12:69611102

94.12 ± 5.56

-0.98 (-1.29 – -0.67)

8.77 × 10-8

RP11-324P9.1

 

2:110449774

84.72 ± 9.25

-1.66 (-2.19 – -1.13)

9.34 × 10-8

BMS1P19

*1 Chromosome number and position of the CpG site.

*2 Gene with evident roles in bone development and/or homeostasis.

*3 Gene with a suggested function in, or relationship with, bone development and/or homeostasis.

DNAm: DNA methylation, SD: standard deviation, CI: confidence interval.

Table 3. CpG sites associated with bone mass (P < 1.0 × 10-7) identified in monocytes

Dataset

Position*1

DNAm level (mean ± SD)

Coefficient

(95% CI)

P value

Gene

Both sexes

2:61108446

0.29 ± 1.38

0.21 (0.15 – 0.26)

2.91 × 10-11

AC010733.4

 

22:47133916

0.84 ± 2.80

0.40 (0.28 – 0.51)

3.73 × 10-10

CERK

 

2:121815071

96.54 ± 4.01

-0.56 (-0.72 – -0.39)

5.95 × 10-10

Y RNA

 

3:188489102

98.19 ± 3.52

-0.46 (-0.60 – -0.32)

3.02 × 10-9

LPP

 

8:8750676

1.11 ± 1.71

0.22 (0.15 – 0.30)

7.51 × 10-9

MFHAS1

 

14:36295936

0.51 ± 1.43

0.19 (0.13 – 0.25)

1.30 × 10-8

RP11-317N8.5/BRMS1L

 

2:195259677

85.29 ± 14.8

-2.22 (-2.89 – -1.55)

2.02 × 10-8

AC018799.1

 

2:69941555

97.83 ± 3.04

-0.39 (-0.52 – -0.26)

2.66 × 10-8

ANXA4

 

3:173900045

85.08 ± 18.96

-2.76 (-3.64 – -1.89)

2.95 × 10-8

NLGN1

 

2:190657692

96.39 ± 3.84

-0.48 (-0.64 – -0.32)

4.33 × 10-8

PMS1

 

12:56474020

0.25 ± 1.11

0.14 (0.09 – 0.18)

6.30 × 10-8

ERBB3*3

 

17:1535968

92.9 ± 13.44

-1.63 (-2.19 – -1.08)

6.43 × 10-8

SCARF1

 

4:99080231

92.82 ± 9.63

-1.41 (-1.86 – -0.95)

6.45 × 10-8

Y RNA

 

8:21530651

94.44 ± 12.01

-1.60 (-2.13 – -1.07)

6.62 × 10-8

GFRA2

 

2:240322590

0.50 ± 1.79

0.22 (0.15 – 0.30)

6.72 × 10-8

HDAC4*3

 

12:129772248

90.49 ± 7.00

-0.85 (-1.14 – -0.56)

8.41 × 10-8

TMEM132D

 

3:16191417

96.63 ± 7.87

-1.00 (-1.34 – -0.66)

9.09 × 10-8

GALNT15

 

14:39639583

1.60 ± 2.61

0.32 (0.21 – 0.43)

9.11 × 10-8

TRAPPC6B

Female

6:35995637

0.37 ± 1.07

0.21 (0.16 – 0.25)

2.44× 10-11

MAPK14*2

 

18:59560980

0.51 ± 1.48

0.27 (0.20 – 0.35)

1.87 × 10-9

RNF152

 

6:163371471

94.89 ± 8.66

-2.42 (-2.96 – -1.87)

2.54 × 10-9

PACRG

 

21:17220599

97.23 ± 4.77

-0.89 (-1.14 – -0.63)

5.16 × 10-9

USP25

 

10:60936039

6.68 ± 5.35

0.98 (0.70 – 1.26)

6.22 × 10-9

PHYHIPL

 

11:497792

96.88 ± 4.27

-0.78 (-1.01 – -0.56)

9.53 × 10-9

RNH1*3

 

4:81105396

0.52 ± 1.96

0.35 (0.25 – 0.45)

1.84 × 10-8

RP11-377G16.2/PRDM8

 

5:65423122

84.35 ± 21.30

-3.75 (-4.89 – -2.60)

3.07 × 10-8

SREK1

 

4:76766071

97.61 ± 3.13

-0.55 (-0.72 – -0.38)

3.34 × 10-8

RP11-556N4.1

 

16:56225078

0.77 ± 2.84

0.50 (0.36 – 0.65)

3.68 × 10-8

RP11-461O7.1

 

20:13765582

0.34 ± 1.48

0.26 (0.18 – 0.34)

3.78 × 10-8

NDUFAF5

 

1:201449877

0.23 ± 0.96

0.17 (0.12 – 0.22)

3.86 × 10-8

CSRP1

 

11:124035292

98.58 ± 3.03

-0.53 (-0.70 – -0.37)

4.21 × 10-8

OR10D1P

 

16:15188395 (cg16603896)

0.96 ± 3.22

0.54 (0.37 – 0.71)

4.60 × 10-8

RP11-72I8.1/PDXDC1

 

2:24660027

92.65 ± 8.78

-1.53 (-2.00 – -1.06)

4.83 × 10-8

NCOA1

 

5:64919917

0.23 ± 0.96

0.16 (0.11 – 0.21)

6.02 × 10-8

TRIM23

 

6:36954391

0.47 ± 1.62

0.29 (0.20 – 0.39)

6.79 × 10-8

MTCH1

 

7:29157080

98.58 ± 3.53

-0.61 (-0.81 – -0.42)

7.31 × 10-8

CPVL

 

2:10338969

97.26 ± 3.76

-0.64 (-0.85 – -0.44)

8.54 × 10-8

C2orf48

 

4:4377945

95.56 ± 3.49

-0.61 (-0.80 – -0.41)

8.58 × 10-8

NSG1

 

1:23857808

0.97 ± 3.41

0.72 (0.50 – 0.93)

9.08 × 10-8

E2F2

 

2:230125683

84.64 ± 9.88

-1.72 (-2.27 – -1.17)

9.61 × 10-8

PID1

Male

10:60540166

87.39 ± 10.86

-2.06 (-2.64 – -1.48)

6.42 × 10-9

BICC1*2

 

1:3463645

93.37 ± 4.07

-0.75 (-0.96 – -0.54)

6.84 × 10-9

MEGF6

 

19:47039387

98.04 ± 4.18

-0.84 (-1.08 – -0.59)

3.90 × 10-8

PPP5D1

 

2:172721808

91.63 ± 11.99

-2.36 (-3.07 – -1.65)

4.68 × 10-8

SLC25A12

 

12:6588863

96.45 ± 3.77

-0.66 (-0.86 – -0.46)

5.54 × 10-8

RP1-102E24.6

 

20:50007014

96.99 ± 4.16

-0.73 (-0.96 – -0.50)

7.50 × 10-8

NFATC2*3

 

14:93054727

92.29 ± 19.21

-3.50 (-4.60 – -2.40)

8.00 × 10-8

RIN3*3

*1 Chromosome number and the position of CpG site.

*2 Gene with evident roles in bone development and/or homeostasis.

*3 Gene with suggested function in, or relationship with, bone development and/or homeostasis. DNAm: DNA methylation, SD: standard deviat