International Journal of Genomics
Volume 2018 (2018), Article ID 6407257, 11 pages
https://doi.org/10.1155/2018/6407257

Research Article

Meta-Analysis of Genome-Wide Association Studies Identifies Novel Functional CpG-SNPs Associated with Bone Mineral Density at Lumbar Spine

Chuan Qiu,¹ Hui Shen,¹ Xiaoying Fu,¹ Chao Xu,¹ and Hongwen Deng1,2

¹Department of Global Biostatistics and Data Science, Center for Bioinformatics and Genomics, School of Public Health and Tropical Medicine, New Orleans 70112, USA
²School of Basic Medical Science, Central South University, Changsha 410013, China

Correspondence should be addressed to Hui Shen; hshen3@tulane.edu

Received 2 May 2018; Accepted 19 July 2018; Published 7 August 2018

Academic Editor: Monika Dmitrzak-Weglarz

Copyright © 2018 Chuan Qiu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Osteoporosis is a serious public health issue, which is mostly characterized by low bone mineral density (BMD). To search for additional genetic susceptibility loci underlying BMD variation, an effective strategy is to focus on testing of specific variants with high potential of functional effects. Single nucleotide polymorphisms (SNPs) that introduce or disrupt CpG dinucleotides (CpG-SNPs) may alter DNA methylation levels and thus represent strong candidate functional variants. Here, we performed a targeted GWAS for 63,627 potential functional CpG-SNPs that may affect DNA methylation in bone-related cells, in five independent cohorts (). By meta-analysis, 9 CpG-SNPs achieved a genome-wide significance level () for association with lumbar spine BMD and additional 15 CpG-SNPs showed suggestive significant () association, of which 2 novel SNPs rs7231498 (NFATC1) and rs7455028 (ESR1) also reached a genome-wide significance level in the joint analysis. Several identified CpG-SNPs were mapped to genes that have not been reported for association with BMD in previous GWAS, such as NEK3 and NFATC1 genes, highlighting the enhanced power of targeted association analysis for identification of novel associations that were missed by traditional GWAS. Interestingly, several genomic regions, such as NEK3 and LRP5 regions, contained multiple significant/suggestive CpG-SNPs for lumbar spine BMD, suggesting that multiple neighboring CpG-SNPs may synergistically mediate the DNA methylation level and gene expression pattern of target genes. Furthermore, functional annotation analyses suggested a strong regulatory potential of the identified BMD-associated CpG-SNPs and a significant enrichment in biological processes associated with protein localization and protein signal transduction. Our results provided novel insights into the genetic basis of BMD variation and highlighted the close connections between genetic and epigenetic mechanisms of complex disease.

1. Introduction

Osteoporosis is a complex disease mainly characterized by low bone mineral density (BMD) and microarchitectural deterioration of bone tissue, which results in an increased risk of bone fragility and susceptibility to fracture [1]. It is an increasingly serious public health issue in the aging population; the prevalence of osteoporosis at lumbar spine in the elderly is over 20% in the United States [2]. Genetic studies have demonstrated that BMD is under strong genetic control, with heritability ranging between 50 and 85% [3, 4]. Genome-wide association studies (GWAS) and meta-analyses of these studies have successfully identified over 250 genetic loci associated with BMDs at different skeletal sites [5–11]. However, these loci explained approximately 12% of BMD variation [11] and the specific functional variants at these loci were generally unknown. To search for additional genetic loci and to enhance our understanding of the biological mechanisms underlying BMD variation, one effective strategy is to focus on testing of specific variants with high potential of functional effects, such as exonic/nonsynonymous variants [5, 12] or variants that may potentially affect regulatory factors [13–16]. Such strategy can alleviate the multiple testing problem of the conventional GWAS approach and consequently enhance the power to identify novel functional variants associated with the phenotype of interest. In addition, because the hypothesis testing is based on SNPs with potential functions, false positive findings may be minimized due to that the information on prior functional evidence is used.

DNA methylation is an essential epigenetic mechanism for the regulation of transcription. It has profound impacts on chromatin structure, genomic imprinting, embryonic development, X-chromosome inactivation, and the pathogenesis of several human genetic disorders [17]. Although epigenetic regulation by DNA methylation is generally thought to be transcriptionally repressed in gene promoters and transcriptionally activated when occurring in gene bodies [18, 19], recent studies suggested a much more complex relationship between DNA methylation and the gene expression pattern. Both positive and negative associations between DNA methylation and gene expression have been revealed across all genomic regions of a gene, and DNA methylation can also modulate alternative RNA splicing via regulation of the RNA Pol II elongation rate [20–24], demonstrating that DNA methylation can have diverse, chromatin cell type- and context-dependent regulatory effects of transcription.

Single nucleotide polymorphisms (SNPs) may introduce or disrupt cytosine-phosphate-guanine dinucleotides (CpG sites), the major substrate for methyl transfer reactions, and therefore dramatically alter the methylation status at the affected loci [25]. These so-called CpG-SNPs have been suggested as an important mechanism through which genetic variants can affect gene function via epigenetics [25, 26]. Shoemaker et al. performed genome-wide allele-specific methylation analysis in 16 human cell lines and found that a significant proportion (38–88%) of allele-specific methylation regions relied on the presence of CpG-SNP variations [27]. Similarly, Zhi et al. conducted genome-wide correlation analysis between genetic variants and DNA methylation levels in human blood CD4+ T cells and found that over 80% of CpG-SNPs were local methylation quantitative trait loci (cis-meQTLs) and CpG-SNPs accounted for over 2/3 of the strongest meQTL signals [28]. The effect of CpG-SNPs often extended beyond the directly affected CpG sites to surrounding regions, likely via correlated proximal methylation patterns and genetic linkage disequilibrium (LD) [25, 28, 29]. These evidences strongly suggested that CpG-SNPs are a crucial type of cis-regulatory polymorphic variants connecting genetic variation to the individual variability in epigenome. By focusing on CpG-SNPs in selected candidate genes, several studies have identified significant associations between CpG-SNPs with human complex disorders, such as breast cancer [30], type 2 diabetes [29], alcohol dependence [31], and suicide attempt in schizophrenia [32], implying that focusing on CpG-SNPs is an efficient strategy to identify novel functional variants underlying human complex disorders/traits.

In this study, we performed a targeted GWAS analysis for BMD on CpG-SNPs. As DNA methylation profiles are often cell-type specific [33], we further narrowed down to CpG-SNPs that are also meQTLs in an osteoclast-lineage cell, specifically, human peripheral blood monocytes (PBMs). PBMs can act as precursors of osteoclasts, produce cytokines important for osteoclast differentiation and function, serve as a major target cell of sex hormones for bone metabolism [34–38], and have been demonstrated as an excellent cell model for studying osteoporosis-related gene/protein expression patterns and their regulatory mechanisms [39–50]. Therefore, our targeted potential functional CpG-SNPs represent prominent candidates that can regulate BMD variation by affecting gene activity via epigenetic mechanisms in bone-related cells.

2. Materials and Methods

2.1. Study Cohorts

The discovery dataset incorporated a total of 5905 subjects from five GWAS, of which three studies were “in-house” studies: (1) Omaha Osteoporosis Study (Caucasian ancestry, ), (2) Kansas City Osteoporosis Study (Caucasian ancestry, ), and (3) China Osteoporosis Study (Han Chinese ancestry, ), and two studies were “external” studies obtained from the Database on Genotypes and Phenotypes (dbGaP): (1) Women’s Health Initiative Observational Study African-American Substudy (African ancestry, ) and (2) Women’s Health Initiative Observational Study Hispanic Substudy (Hispanic ancestry, ). The basic characteristics of the five study cohorts were shown in Supplementary Table 1. All studies were reviewed and approved from respective institutional review boards, and each eligible participant provided written informed consent for enrolment. The replication dataset included the summary statistics for the association of approximately 10 million SNPs with BMD by the Genetic Factors for Osteoporosis Consortium (GEFOS) [5, 8]. To our knowledge, it is the largest GWAS meta-analysis dataset for BMD association to date in the bone field [5, 8].

2.2. Selecting Potential Functional CpG-SNPs

The CpG-SNPs that are potentially functional in PBMs were selected according to the following steps:

(1)We identified CpG-SNPs in the human genome by interrogating the extensive catalog of common and rare genetic variants from the 1000 Genomes reference panel [51] and our in-house whole-genome high-coverage deep resequencing study [52]. A SNP was defined as a CpG-SNP if it introduces or disrupts a CpG site. A total of 3,363,517 CpG-SNPs was identified throughout the human genome.(2)We retrieved 39,859 PBM meQTLs at a stringent significance threshold (FDR < 0.001) from the previous study that assessed the association of over 7 million SNPs with methylome of PBMs in 200 unrelated individuals [53]. We then used SNiPA [54] to identify proxy SNPs in strong LD with retrieved PBM meQTLs. The search was depended on genotype information from the 1000 Genomes Project with the European samples [51]. The inclusion criteria for proxy SNPs were set as a pairwise ² threshold > 0.9 and a distance limit of 10 kb from the query meQTL. A total of 175,710 potential PBM DNA methylation-associated SNPs (reported meQTLs and proxy of meQTLs) were identified.(3)By finding common SNPs between the CpG-SNPs and PBM DNA methylation-associated SNPs, we identified a total of 68,041 CpG-SNPs that are potentially functional in PBMs.

2.3. BMD Measurements

The lumbar spine BMD was determined by either the Hologic Inc. (Bedford, MA, USA) or GE Lunar Corp. (Madison, WI, USA) dual-energy X-ray absorptiometry (DXA) scanner following the respective manufacturer’s scan protocols. For each GWAS, multiple potential covariates such as scanner ID, sex, height, weight, age, and age² were screened using a forward stepwise linear regression. The significant covariates were used to adjust for raw BMD measurements. Correction of potential population stratification was performed with principal component analysis (PCA), and the top five PCs (i.e., PC1–PC5) were also included as covariates. Residual scores of adjusted phenotypes were normalized by inverse quantile of the standard normal distribution, which was analyzed subsequently.

2.4. Genotyping and Quality Control

For each GWAS, genome-wide genotyping was performed by either Affymetrix Inc. (Santa Clara, CA, USA) or Illumina Inc. (San Diego, CA, USA) high-density SNP genotyping platforms following respective manufacturer’s assay protocols. Quality control was implemented by PLINK (http://pngu.mgh.harvard.edu/~purcell/plink/) with the following criteria: individual missingness < 5%, SNP with successful call rate > 95%, and Hardy-Weinberg equilibrium value > 1.0 × 10⁻⁵. PCs derived from genome-wide genotyping analysis were used to monitor the population outliers.

2.5. Genotype Imputation

To allow for the merging of datasets from different types of genotyping platform to obtain higher depth of genome coverage, we performed extensive genotype imputation analysis. Generally, haplotype inference of each GWA study was initially phased by a Markov Chain Haplotyping algorithm (MACH) [55] and Minimac [56] was then used to impute genotypes at untyped variants based on haplotype data from the 1000 Genomes reference panel [51]. For each GWA study, the haplotype reference panel of relevant population was used to impute genotypes at untyped variants. SNPs with imputation quality score (²) > 0.3 and minor allele frequency (MAF) > 0.05 in no less than 2 studies were retained in the subsequent analyses. Imputation with the 1000 Genomes Project reference panels generated genotype data for more than 11.2 million SNPs. Among the 68,041 potential functional CpG-SNPs, 63,627 CpG-SNPs had qualified genotype data (genotyped + imputed) and thus were tested in the following GWAS meta-analyses.

2.6. Association Tests and Meta-Analyses

For each GWAS, we test the association between directly typed/imputed SNPs and lumbar spine BMD using an additive genetic model. The association of unrelated subjects in each GWAS was tested by fitting a linear regression model with MACH2QTL [55] in which allele dosage was considered as a phenotype predictor. The genomic inflation factor (λGC) [57] was also estimated for each individual GWAS. We performed meta-analysis using software METAL [58] which based on weights proportional to the square root of the number of subjects in each sample, and between-study heterogeneity was estimated by Cochran’s statistic and ². Genome-wide significance threshold was defined as a value < 7.86 × 10⁻⁷ (Bonferroni correction for testing 63,627 selected CpG-SNPs).

2.7. Function Annotation of the CpG-SNPs

CpG-SNPs were annotated with SNPnexus [59] based on reference genome GRCh37 and assigned to candidate genes (±2 kb upstream and downstream). In order to test the potential functional importance of the identified CpG-SNPs, we applied HaploReg [60] to annotate selected CpG-SNPs to enhancer histone marks (H3K4me1/H3K27ac) across diverse tissue/cell types from the Roadmap Epigenomics Projects and test the effect of SNPs on changing the regulatory motifs and the effect of SNPs on the regulation of gene expression of target genes. We employed the software GOEAST [61] to identify significant gene ontology terms among genes associated with identified novel functional CpG-SNPs in lumbar spine.

3. Results

In this study, we identified 68,041 potential functional CpG-SNPs that may both affect DNA methylation by introducing or disrupting CpG sites and influence DNA methylation levels in human PBMs. Interestingly, although over 50% of these potential functional CpG-SNPs were mapped to introns, we observed a significant enrichment of potential functional CpG-SNPs in 5/3-UTR regions (fold change > 2) and underrepresentation in intergenic regions (Supplementary Figure 1), when comparing to the overall profile of CpG-SNPs in the human genome.

We successfully obtained genotype data for 63,627 potential functional CpG-SNPs and carried out targeted association studies in five independent GWAS cohorts with a total of 5905 subjects. The estimates of genomic inflation factor λGC ranged from 0.97 to 1.02 in individual GWAS. By performing meta-analysis combining the five GWAS datasets, we identified 9 CpG-SNPs that were significantly associated with lumbar spine BMD at a genome-wide significance level (), including 5 novel SNPs rs689179 ( value = 2.68 × 10⁻⁷), rs576118 ( value = 2.70 × 10⁻⁷), rs471966 ( value = 3.29 × 10⁻⁷), rs640569 ( value = 4.04 × 10⁻⁷), and rs667126 ( value = 7.80 × 10⁻⁷) in LRP5 gene and one SNP rs9535889 in novel gene NEK3 ( value = 7.55 × 10⁻⁷). We also confirmed 3 previously reported loci (rs525592, rs1784235, and rs497261) in LRP5 gene (Figure 1 and Table 1). In addition, 15 CpG-SNPs achieved a suggestive significance level () for association with lumbar spine BMD (Figure 1 and Table 1). We then performed in silico replication for the identified 24 significant/suggestive CpG-SNPs in the GEFOS cohort [5, 8] and successfully replicated ( value < 0.05) 14 CpG-SNPs (Table 1). Subsequently, a joint analysis of both the discovery and replication studies identified 2 additional novel CpG-SNPs associated with lumbar spine BMD at a genome-wide significance level (Table 1) including the SNP rs7455028 ( value = 1.18 × 10⁻⁷) in ESR1 gene and the SNP rs7231498 ( value = 7.18 × 10⁻⁷) in NFATC1 gene. A number of the significant/suggestive CpG-SNPs were clustered into the genomic regions encompassing NEK3 and LRP5 genes (Figure 2 and Supplementary Figure 2). These clustered CpG-SNPs are in high LD and therefore, may represent the same functional loci that synergistically mediate the DNA methylation and/or gene expression of their target genes.

Figure 1: Circular Manhattan plot picturing the −log₁₀ ( values) of meta-analysis results for lumbar spine BMD. CpG-SNPs were plotted according to the chromosomal location. The blue and red circular lines indicate the threshold for suggestive significant ( value = 5.00 × 10⁻⁵) and significant SNPs ( value = 7.86 × 10⁻⁷), respectively.

Table 1: Significant/suggestive CpG-SNPs for lumbar spine BMD.

CpG-SNP
Chr
Position
Alleles
Nearest gene
Feature
Meta value
GEFOS value
Joint value

rs2941741
6
152008982
G/A
ESR1
Intronic
6.50E − 06
1.21E − 08
2.45E − 12

rs3020333
6
152010254
A/G
ESR1
Intronic, 5 upstream
7.17E − 06
1.57E − 09
3.73E − 13

rs7455028
6
152034386
C/T
ESR1
Intronic
4.57E − 05
0.00013
1.18E − 07

rs13254554
8
120010805
T/C
COLEC10/TNFRSF11B
Intronic
8.99E − 07
1.10E − 19
5.79E − 24

rs2220189
8
120007708
C/G
COLEC10/TNFRSF11B
Intronic
1.53E − 06
4.25E − 20
3.84E − 24

rs525592
11
68195104
C/T
LRP5
Intronic
1.86E − 07
8.69E − 11
6.41E − 16

rs689179
11
68179166
A/G
LRP5
Intronic
2.68E − 07
NA
NA

rs576118
11
68177708
G/A
LRP5
Intronic
2.70E − 07
3.81E − 06
2.95E − 11

rs471966
11
68173861
C/T
LRP5
Intronic
3.29E − 07
NA
NA

rs1784235
11
68185500
C/T
LRP5
Intronic
3.92E − 07
2.95E − 08
3.82E − 13

rs640569
11
68184820
A/G
LRP5
Intronic
4.04E − 07
2.17E − 08
2.92E − 13

rs667126
11
68177728
C/T
LRP5
Intronic
7.80E − 07
NA
NA

rs497261
11
68192244
T/C
LRP5
Intronic
7.83E − 07
1.86E − 11
5.79E − 16

rs314751
11
68179560
C/T
LRP5
Intronic
1.31E − 06
1.20E − 11
6.24E − 16

rs23691
11
68178668
G/A
LRP5
Intronic
1.33E − 06
1.17E − 11
6.18E − 16

rs531163
11
68194496
A/G
LRP5
Intronic
1.35E − 06
9.05E − 11
4.60E − 15

rs9535889
13
52733634
C/G
NEK3
Intronic, 5 upstream, 5 UTR
7.55E − 07
0.469098
5.61E − 06

rs3783242
13
52717950
C/T
NEK3
Intronic, 3 downstream
2.41E − 06
NA
NA

rs9526841
13
52726476
A/G
NEK3
Intronic, 3 downstream
6.64E − 06
0.325285
3.03E − 05

rs2897976
13
52715944
G/A
NEK3
Intronic
9.71E − 06
0.471396
6.08E − 05

rs9526843
13
52730056
C/T
NEK3
Intronic
1.17E − 05
NA
NA

rs2408609
13
52714043
C/T
NEK3
Intronic
1.32E − 05
0.546887
9.26E − 05

rs2408611
13
52709742
G/A
NEK3
Intronic, 5 upstream
2.50E − 05
NA
NA

rs7231498
18
77189387
A/G
NFATC1
Intronic
4.22E − 05
0.000944
7.18E − 07

Note: CpG-SNPs reached a genome-wide significance level ( value ≤ 7.86 × 10⁻⁷) in discovery meta-analysis and/or joint analysis of discovery, and replication studies are marked in bold. Gene/CpG-SNP reported in previous GWAS for BMD is marked in italics. NA: SNPs were not available in the GEFOS 2015 data release. This result was based on the GEFOS 2012 data release because this SNP is not available in the 2015 release.

Figure 2: A regional association plots of significant/suggestive CpG-SNPs at NEK3 regions. Genes and expressed sequence tags (ESTs) within the region are shown in the lower panel, and the unbroken blue line indicates the recombination rate within the region. Each filled circle represents the value for one SNP in the meta-analysis, with the top SNP rs9535889 shown in purple and SNPs in the region colored depending on their degree of LD (²) with rs9535889. LD was estimated by LocusZoom [80] on the basis of CEU (Utah residents of Northern and Western European ancestry) HapMap haplotype data.

To further explore the potential functional significance of the identified significant/suggestive CpG-SNPs, we annotated these CpG-SNPs to various chromatin states and other possible regulatory elements with data from Roadmap Epigenomics and GTEx projects through the HaploReg program [60]. The chromatin state and histone modification data suggested the evidence of regulatory potential in the identified CpG-SNPs. 20 CpG-SNPs altered the regulatory motif, along with 14 CpG-SNPs involved enhancer histone markers. Notably, the novel SNPs rs9535889, rs9526841, and rs2408611 in NEK3 gene were all located in regions with strong transcription and enhancer activities in PBMs as well as various other tissues and cell types (Table 2), highlighting strong regulatory potential of these CpG-SNPs. In addition, many identified CpG-SNPs may affect binding of various transcription factors and have numerous reported eQTL evidences in various tissue/cell types (Table 2). We also conducted gene ontology analysis for the genes related to the identified CpG-SNPs and revealed significant enrichment of biological processes which are closely associated to protein localization and protein signal transduction (Table 3), such as protein localization to plasma membrane/cell periphery and regulation of Ras/Rho protein signal transduction gene ontology terms.

Table 2: Functional annotation of significant/suggestive CpG-SNPs.

CpG-SNPs
Nearest gene
Chromatin state in PBMs¹
Tissues/cells with enhancer histone marks (H3K4me1/H3K27ac)
Motifs changed
eQTL hits

rs2941741
ESR1
Quiescent/low

rs3020333
ESR1
Quiescent/low
Liver
Pou2f2

rs7455028
ESR1
Quiescent/low

5 altered motifs

rs13254554
COLEC10/TNFRSF11B
Quiescent/low

TCF12, p53
3 hits

rs2220189
COLEC10/TNFRSF11B
Quiescent/low
7 tissues

2 hits

rs525592
LRP5
Quiescent/low

4 altered motifs
2 hits

rs689179
LRP5
Quiescent/low
6 tissues
5 altered motifs
1 hit

rs576118
LRP5
Quiescent/low
IPSC, muscle, heart
TAL1
1 hit

rs471966
LRP5
Quiescent/low
8 tissues
8 altered motifs
2 hits

rs1784235
LRP5
Quiescent/low
Blood
AP-2, ELF1, Rad21
1 hit

rs640569
LRP5
Quiescent/low
Blood
Irf, Pax-4, Pou2f2
1 hit

rs667126
LRP5
Quiescent/low
IPSC, muscle, heart
9 altered motifs
2 hits

rs497261
LRP5
Quiescent/low
Muscle
Pax-5, Smad
5 hits

rs314751
LRP5
Quiescent/low
6 tissues
4 altered motifs
4 hits

rs23691
LRP5
Quiescent/low
6 tissues
5 altered motifs
4 hits

rs531163
LRP5
Quiescent/low

7 altered motifs
3 hits

rs9535889
NEK3
Active TSS
24 tissues²
Rad21, SP1, TATA
72 hits

rs3783242
NEK3
Quiescent/low

CDP, Pou2f2
72 hits

rs9526841
NEK3
Strong transcription

HIF1, RFX5, TCF11::MafG
62 hits

rs2897976
NEK3
Quiescent/low

79 hits

rs9526843
NEK3
Quiescent/low
Intestine

73 hits

rs2408609
NEK3
Quiescent/low

6 altered motifs
78 hits

rs2408611
NEK3
Strong transcription

4 altered motifs
80 hits

rs7231498
NFATC1
Weak transcription
Blood
7 altered motifs

Note: ¹Chromatin state information was retrieved using a 15-state model from the Roadmap Epigenomics Project based on the 5 core histone marks. ²Tissues/cells with promoter histone marks (H3K4me3/H3K9ac). Abbreviation: IPSC: induced pluripotent stem cells.

Table 3: The top ten most significant GO terms enriched for BMD-associated CpG-SNPs.

GOID
Term
Log odds ratio
value

GO:0072659
Protein localization to plasma membrane
5.42
4.06E − 16

GO:1990778
Protein localization to cell periphery
5.42
4.06E − 16

GO:0007009
Plasma membrane organization
4.96
3.08E − 14

GO:0035023
Regulation of Rho protein signal transduction
4.57
1.21E − 12

GO:0046578
Regulation of Ras protein signal transduction
4.01
1.96E − 10

GO:0072657
Protein localization to membrane
3.94
3.45E − 10

GO:0010256
Endomembrane system organization
3.93
3.75E − 10

GO:0000904
Cell morphogenesis involved in differentiation
2.92
4.50E − 08

GO:0051056
Regulation of small GTPase-mediated signal transduction
3.29
7.92E − 08

GO:0008295
Spermidine biosynthetic process
7.99
2.55E − 07

Note: GO enrichment analysis was performed in candidate genes annotated to BMD-associated CpG-SNPs ( value < 1.0 × 10⁻⁴).

4. Discussion

Our study represents the first targeted GWAS testing CpG-SNPs that are potentially functional in bone-related cells for association with BMD variation. As epigenomic and transcriptomic profiles are often tissue-/cell-type specific, we speculated that only a subset of CpG-SNPs in the human genome will have functional impact on DNA methylation levels in specific tissues/cells. Therefore, it is necessary to select out those CpG-SNPs that are potentially functional in disease-/trait-related tissues/cells when performing CpG-SNP-focused association studies. One reasonable and efficient filtering strategy is to leverage the enormous available meQTL data in diverse tissues/cells. Unfortunately, meQTL data in skeletal cells were scarce; therefore, we used meQTL data from PBMs to select out 68,041 candidate CpG-SNPs that may be functional in regulating bone mass, considering the direct and close connections between PBMs and bone metabolism. These potential functional CpG-SNPs were enriched in 5/3-UTR regions but underrepresented in intergenic regions (Supplementary Figure 1). This result is largely in line with the recent findings in tissue-/cell-type specific DNA methylation profiles, suggesting that methylation-mediated regulatory effects often occur beyond the promoter areas [18, 62].

By using the data from five independent GWAS cohorts and the summary statistics from the GEFOS study, we identified significant/suggestive associations for 24 CpG-SNPs with lumbar spine BMD. These BMD-associated CpG-SNPs were mapped to six genes; some of which have not been reported for association with BMD in previous GWAS, such as NEK3 and NFATC1 genes. Our finding highlighted the enhanced power of targeted association analysis for identification of novel associations that were missed by traditional GWAS. Interestingly, several genomic regions, such as LRP5 and NEK3 regions, contained multiple significant/suggestive CpG-SNPs, suggesting that multiple neighboring CpG-SNPs may synergistically mediate the DNA methylation and gene expression of the target genes. This is consistent with the fact that methylation signals among neighboring CpG sites are often strongly correlated and regulatory elements that were mediated by methylation usually extend across various genomic regions [63]. LRP5 gene encodes a transmembrane protein which acts as a receptor for low-density lipoprotein. This transmembrane receptor initializes the process of receptor-mediated endocytosis by binding and internalizing their corresponding ligands [64]. It is well known for the critical role in bone homeostasis and several skeletal disorders [65]. Several common genetic variants of LRP5 gene have been demonstrated as potential risk factors in osteoporosis and fracture by previous GWAS [66, 67]. For example, gain of functional variations in LRP5 gene leads to extremely high BMD [64] and loss of functional variations in LRP5 gene results in osteoporosis-pseudoglioma syndrome [68]. Interestingly, the recent study showed that the differentiation of monocytes can be negatively regulated by LRP5 gene through abrogation of the Wnt pathway which has an essential role in bone remodeling in both physiological and pathological conditions [69]. The other interesting gene is the NEK3, which also contained several significant/suggestive CpG-SNPs and enriched with strong transcription and enhancer histone modification marks in PBMs and a variety of other tissue/cell types. NEK3 gene encodes a member of the NimA-related serine/threonine kinases [70]. These kinases have been implicated as the significant regulators of cell migration [71] and also regulate microtubule acetylation in neurons [72]. Although most of these CpG-SNPs were annotated to introns of the NEK3 gene, the eQTL data from GTEx project suggested that these CpG-SNPs were strongly associated with the expression of NEK3 gene in diverse tissues. Notably, the previous study [73] that assessed the association of over 675,000 SNPs with transcriptome of PBMs in 1490 unrelated individuals showed that SNP rs2408611 in NEK3 gene has a strong cis-eQTL effect in PBM. This evidence may support that CpG-SNP-mediated epigenomic alterations may be an important mechanism underlying the association between NEK3 and BMD variation. However, its function in other tissues, including bone, remains largely uncharacterized. Another interesting gene is NFATC1. This gene encodes a transcription factor involved in T cell maturation. Importantly, NFATC1 can also regulate activity of a number of osteoclast-specific enzymes and/or other molecules, such as osteoclast-associated receptor, TRAP, calcitonin receptor, and cathepsin K through cooperation with MITF and c-Fos [74–77]. The important role of this gene in differentiation of osteoclast has been well established by several studies performed on genetically modified mutant mice [78, 79]. For example, Winslow et al. [78] identified that the transgenic mice generated by crossing NFATC1-knockout mice with mice that express Tie2 promoter-driven NFATC1 exhibit an osteopetrotic bone phenotype, which may result from a severe defect in the osteoclastogenesis process. Therefore, understanding the molecular basis underlying the functional regulation of NFATC1 in osteoclasts may provide novel therapeutic strategies for bone diseases.

Several potential limitations of this study should be concerned and addressed in the future. First, the selection of an appropriate cell model is crucial. Due to the limited meQTL studies in the bone cell models, here, we focused on CpG-SNPs that are also meQTLs in an osteoclast-lineage cell, specifically, human PBMs. Although PBMs act as precursors of osteoclasts and act as the major target cells of sex hormones for bone metabolism, the ideal model cells for the osteoporosis study are bone cells, such as osteoblast, osteoclast, and osteocyte. Second, the results of functional annotation exclusively depend on computationally predicted regulation features and further experimental validation should be conducted to confirm the biological significance of these potential functional CpG-SNPs.

In summary, we performed a targeted GWAS analysis for potential functional CpG-SNPs and identified 2 novel BMD-associated genes, NEK3 and NFATC1. Our results highlighted the power of targeted analysis of potential functional variants for the identification of novel disease susceptibility loci that have been missed by a conventional GWAS approach. More importantly, our findings suggested that CpG-SNP-mediated DNA methylation changes may be a crucial biological mechanism to be considered in the interpretation of associations between common genetic variants, epigenetic process, and phenotypes of human diseases.

Data Availability

The chromatin data used to support the findings of this study have been deposited in the Roadmap Epigenomics Project repository. The genotype data used to support the findings of this study are available from the corresponding author upon request.

Disclosure

This manuscript was not prepared in collaboration with investigators of the WHI, has not been reviewed and/or approved by WHI investigators, and does not necessarily reflect the opinions of the WHI investigators or the NHLBI.

Conflicts of Interest

The authors declare that there is no conflict of interest.

Authors’ Contributions

Chuan Qiu and Hui Shen contributed equally to the work.

Acknowledgments

This study was partially supported or benefited by grants from the National Institutes of Health (R01AR059781, P20GM109036, R01MH107354, R01MH104680, R01GM109068, R01AR069055, and U19AG055373), the Franklin D. Dickson/Missouri Endowment, the Edward G. Schlieder Endowment, and the Drs. W. C. Tsai and P. T. Kung Professorship in Biostatistics from Tulane University. The Women’s Health Initiative (WHI) program is funded by the National Heart, Lung, and Blood Institute, National Institutes of Health, U.S. Department of Health and Human Services through Contracts N01WH22110, 24152, 32100-2, 32105-6, 32108-9, 32111-13, 32115, 32118-32119, 32122, 42107-26, 42129-32, and 44221. WHI Population Architecture Using Genomics and Epidemiology (PAGE) is funded through the NHGRI Population Architecture Using Genomics and Epidemiology (PAGE) network (Grant no. U01 HG004790). Assistance with phenotype harmonization, SNP selection, data cleaning, meta-analyses, data management and dissemination, and general study coordination was provided by the PAGE Coordinating Center (U01HG004801-01). The authors thank the authors who generously shared their data.

Supplementary Materials

Supplementary Table 1: basic characteristics of the studied samples. Supplementary Figure 1: distribution of CpG-SNPs in distinct genomic features. Supplementary Figure 2: a regional association plots of significant/suggestive CpG-SNPs at LRP5 regions. (Supplementary Materials)

References

J. A. Kanis, “Diagnosis of osteoporosis and assessment of fracture risk,” The Lancet, vol. 359, no. 9321, pp. 1929–1936, 2002.

N. C. Wright, A. C. Looker, K. G. Saag et al., “The recent prevalence of osteoporosis and low bone mass in the United States based on bone mineral density at the femoral neck or lumbar spine,” Journal of Bone and Mineral Research, vol. 29, no. 11, pp. 2520–2526, 2014.

S. H. Ralston and B. de Crombrugghe, “Genetic regulation of bone mass and susceptibility to osteoporosis,” Genes & Development, vol. 20, no. 18, pp. 2492–2506, 2006.

S. H. Ralston and A. G. Uitterlinden, “Genetics of osteoporosis,” Endocrine Reviews, vol. 31, no. 5, pp. 629–662, 2010.

H. F. Zheng, V. Forgetta, Y. H. Hsu et al., “Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture,” Nature, vol. 526, no. 7571, pp. 112–117, 2015.

L. Zhang, H. J. Choi, K. Estrada et al., “Multistage genome-wide association meta-analyses identified two new loci for bone mineral density,” Human Molecular Genetics, vol. 23, no. 7, pp. 1923–1933, 2014.

D. L. Koller, H. F. Zheng, D. Karasik et al., “Meta-analysis of genome-wide studies identifies WNT16 and ESR1 SNPs associated with bone mineral density in premenopausal women,” Journal of Bone and Mineral Research, vol. 28, no. 3, pp. 547–558, 2013.

K. Estrada, U. Styrkarsdottir, E. Evangelou et al., “Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture,” Nature Genetics, vol. 44, no. 5, pp. 491–501, 2012.

Shen, C. Qiu, and H.-W. Deng, “Genetics of osteoporotic fracture,” Orthopedic Research and Reviews, vol. 3, pp. 11–21, 2011.

Y. Guo, L. J. Tan, S. F. Lei et al., “Genome-wide association study identifies ALDH7A1 as a novel susceptibility gene for osteoporosis,” PLoS Genetics, vol. 6, no. 1, article e1000806, 2010.

J. P. Kemp, J. A. Morris, C. Medina-Gomez et al., “Identification of 153 new loci associated with heel bone mineral density and functional involvement of GPC6 in osteoporosis,” Nature Genetics, vol. 49, no. 10, pp. 1468–1475, 2017.

H. R. Kim, N. M. Duc, and K. Y. Chung, “Comprehensive analysis of non-synonymous natural variants of G protein-coupled receptors,” Biomolecules & Therapeutics, vol. 26, no. 2, pp. 101–108, 2018.

T. Niu, N. Liu, X. Yu et al., “Identification of IDUA and WNT16 phosphorylation-related non-synonymous polymorphisms for bone mineral density in meta-analyses of genome-wide association studies,” Journal of Bone and Mineral Research, vol. 31, no. 2, pp. 358–368, 2016.

S. F. Lei, C. J. Papasian, and H. W. Deng, “Polymorphisms in predicted miRNA binding sites and osteoporosis,” Journal of Bone and Mineral Research, vol. 26, no. 1, pp. 72–78, 2011.

O. I. Koues, R. A. Kowalewski, L. W. Chang et al., “Enhancer sequence variants and transcription-factor deregulation synergize to construct pathogenic regulatory circuits in B-cell lymphoma,” Immunity, vol. 42, no. 1, pp. 186–198, 2015.

A. Takata, I. Ionita-Laza, J. A. Gogos, B. Xu, and M. Karayiorgou, “De novo synonymous mutations in regulatory elements contribute to the genetic etiology of autism and schizophrenia,” Neuron, vol. 89, no. 5, pp. 940–947, 2016.

M. Weber, I. Hellmann, M. B. Stadler et al., “Distribution, silencing potential and evolutionary impact of promoter DNA methylation in the human genome,” Nature Genetics, vol. 39, no. 4, pp. 457–466, 2007.

P. A. Jones, “Functions of DNA methylation: islands, start sites, gene bodies and beyond,” Nature Reviews Genetics, vol. 13, no. 7, pp. 484–492, 2012.

M. P. Ball, J. B. Li, Y. Gao et al., “Targeted and genome-scale strategies reveal gene-body methylation signatures in human cells,” Nature Biotechnology, vol. 27, no. 4, pp. 361–368, 2009.

S. Chandra, C. Baribault, M. Lacey, and M. Ehrlich, “Myogenic differential methylation: diverse associations with chromatin structure,” Biology, vol. 3, no. 2, pp. 426–451, 2014.

M. Ehrlich and M. Lacey, “DNA methylation and differentiation: silencing, upregulation and modulation of gene expression,” Epigenomics, vol. 5, no. 5, pp. 553–568, 2013.

Y. Liu, J. Ding, L. M. Reynolds et al., “Methylomics of gene expression in human monocytes,” Human Molecular Genetics, vol. 22, no. 24, pp. 5065–5074, 2013.

J. Wan, V. F. Oliver, G. Wang et al., “Characterization of tissue-specific differential DNA methylation suggests distinct modes of positive and negative gene expression regulation,” BMC Genomics, vol. 16, no. 1, p. 49, 2015.

M. D. Long, D. J. Smiraglia, and M. J. Campbell, “The genomic impact of DNA CpG methylation on gene expression; relationships in prostate cancer,” Biomolecules, vol. 7, no. 4, 2017.

J. Gertz, K. E. Varley, T. E. Reddy et al., “Analysis of DNA methylation in a three-generation family reveals widespread genetic influence on epigenetic regulation,” PLoS Genetics, vol. 7, no. 8, article e1002228, 2011.

J. T. Bell, A. A. Pai, J. K. Pickrell et al., “DNA methylation patterns associate with genetic and gene expression variation in HapMap cell lines,” Genome Biology, vol. 12, no. 1, article R10, 2011.

R. Shoemaker, J. Deng, W. Wang, and K. Zhang, “Allele-specific methylation is prevalent and is contributed by CpG-SNPs in the human genome,” Genome Research, vol. 20, no. 7, pp. 883–889, 2010.

D. Zhi, S. Aslibekyan, M. R. Irvin et al., “SNPs located at CpG sites modulate genome-epigenome interaction,” Epigenetics, vol. 8, no. 8, pp. 802–806, 2013.

T. A. Dayeh, A. H. Olsson, P. Volkov, P. Almgren, T. Ronn, and C. Ling, “Identification of CpG-SNPs associated with type 2 diabetes and differential DNA methylation in human pancreatic islets,” Diabetologia, vol. 56, no. 5, pp. 1036–1046, 2013.

S. Harlid, M. I. L. Ivarsson, S. Butt et al., “A candidate CpG SNP approach identifies a breast cancer associated ESR1-SNP,” International Journal of Cancer, vol. 129, no. 7, pp. 1689–1698, 2011.

M. M. Taqi, I. Bazov, H. Watanabe et al., “Prodynorphin CpG-SNPs associated with alcohol dependence: elevated methylation in the brain of human alcoholics,” Addiction Biology, vol. 16, no. 3, pp. 499–509, 2011.

A. Bani-Fatemi, V. F. Goncalves, C. Zai et al., “Analysis of CpG SNPs in 34 genes: association test with suicide attempt in schizophrenia,” Schizophrenia Research, vol. 147, no. 2-3, pp. 262–268, 2013.

K. E. Varley, J. Gertz, K. M. Bowling et al., “Dynamic DNA methylation across diverse human cell lines and tissues,” Genome Research, vol. 23, no. 3, pp. 555–567, 2013.

Y. Fujikawa, J. M. Quinn, A. Sabokbar, J. O. McGee, and N. A. Athanasou, “The human osteoclast precursor circulates in the monocyte fraction,” Endocrinology, vol. 137, no. 9, pp. 4058–4060, 1996.

S. Higuchi, N. Tabata, M. Tajima et al., “Induction of human osteoclast-like cells by treatment of blood monocytes with anti-fusion regulatory protein-1/CD98 monoclonal antibodies,” Journal of Bone and Mineral Research, vol. 13, no. 1, pp. 44–49, 1998.

A. Matayoshi, C. Brown, J. F. DiPersio et al., “Human blood-mobilized hematopoietic precursors differentiate into osteoclasts in the absence of stromal cells,” Proceedings of the National Academy of Sciences of the United States of America, vol. 93, no. 20, pp. 10785–10790, 1996.

Y. Komano, T. Nanki, K. Hayashida, K. Taniguchi, and N. Miyasaka, “Identification of a human peripheral blood monocyte subset that differentiates into osteoclasts,” Arthritis Research & Therapy, vol. 8, no. 5, article R152, 2006.

R. Lari, P. D. Kitchener, and J. A. Hamilton, “The proliferative human monocyte subpopulation contains osteoclast precursors,” Arthritis Research & Therapy, vol. 11, no. 1, article R23, 2009.

Y. Mori, S. Tsuji, M. Inui et al., “Inhibitory immunoglobulin-like receptors LILRB and PIR-B negatively regulate osteoclast development,” The Journal of Immunology, vol. 181, no. 7, pp. 4742–4751, 2008.

C. S. Soltanoff, S. Yang, W. Chen, and Y. P. Li, “Signaling networks that control the lineage commitment and differentiation of bone cells,” Critical Reviews™ in Eukaryotic Gene Expression, vol. 19, no. 1, pp. 1–46, 2009.

B. Sung, A. Murakami, B. O. Oyajobi, and B. B. Aggarwal, “Zerumbone abolishes RANKL-induced NF-κB activation, inhibits osteoclastogenesis, and suppresses human breast cancer-induced bone loss in athymic nude mice,” Cancer Research, vol. 69, no. 4, pp. 1477–1484, 2009.

M. Durand, S. V. Komarova, A. Bhargava et al., “Monocytes from patients with osteoarthritis display increased osteoclastogenesis and bone resorption: the in vitro osteoclast differentiation in arthritis study,” Arthritis and Rheumatism, vol. 65, no. 1, pp. 148–158, 2013.

R. Leung, K. Cuddy, Y. Wang, J. Rommens, and M. Glogauer, “Sbds is required for Rac2-mediated monocyte migration and signaling downstream of RANK during osteoclastogenesis,” Blood, vol. 117, no. 6, pp. 2044–2053, 2011.

R. Matsubara, T. Kukita, Y. Ichigi et al., “Characterization and identification of subpopulations of mononuclear preosteoclasts induced by TNF-α in combination with TGF-β in rats,” PLoS One, vol. 7, no. 10, article e47930, 2012.

K. Y. Park, W. A. Li, and M. O. Platt, “Patient specific proteolytic activity of monocyte-derived macrophages and osteoclasts predicted with temporal kinase activation states during differentiation,” Integrative Biology, vol. 4, no. 12, pp. 1459–1469, 2012.

L. J. Guo, L. Liao, L. Yang, Y. Li, and T. J. Jiang, “MiR-125a TNF receptor-associated factor 6 to inhibit osteoclastogenesis,” Experimental Cell Research, vol. 321, no. 2, pp. 142–152, 2014.

P. T. Yang, X. H. Meng, Y. Yang, and W. G. Xiao, “Inhibition of osteoclast differentiation and matrix metalloproteinase production by CD4+CD25+ T cells in mice,” Osteoporosis International, vol. 24, no. 3, pp. 1113-1114, 2013.

C. R. Farber, “Identification of a gene module associated with BMD through the integration of network analysis and genome-wide association data,” Journal of Bone and Mineral Research, vol. 25, no. 11, pp. 2359–2367, 2010.

M. Kotani, J. Kikuta, F. Klauschen et al., “Systemic circulation and bone recruitment of osteoclast precursors tracked by using fluorescent imaging techniques,” Journal of Immunology, vol. 190, no. 2, pp. 605–612, 2013.

Y. Zhou, H. W. Deng, and H. Shen, “Circulating monocytes: an appropriate model for bone-related study,” Osteoporosis International, vol. 26, no. 11, pp. 2561–2572, 2015.

R. M. Durbin, D. L. Altshuler, R. M. Durbin et al., “A map of human genome variation from population-scale sequencing,” Nature, vol. 467, no. 7319, pp. 1061–1073, 2010.

H. Shen, J. Li, J. Zhang et al., “Comprehensive characterization of human genome variation by high coverage whole-genome sequencing of forty four Caucasians,” PLoS One, vol. 8, no. 4, article e59494, 2013.

L. Chen, B. Ge, F. P. Casale et al., “Genetic drivers of epigenetic and transcriptional variation in human immune cells,” Cell, vol. 167, no. 5, pp. 1398–1414.e24, 2016.

M. Arnold, J. Raffler, A. Pfeufer, K. Suhre, and G. Kastenmuller, “SNiPA: an interactive, genetic variant-centered annotation browser,” Bioinformatics, vol. 31, no. 8, pp. 1334–1336, 2015.

Y. Li, C. J. Willer, J. Ding, P. Scheet, and G. R. Abecasis, “MaCH: using sequence and genotype data to estimate haplotypes and unobserved genotypes,” Genetic Epidemiology, vol. 34, no. 8, pp. 816–834, 2010.

C. Fuchsberger, G. R. Abecasis, and D. A. Hinds, “minimac2: faster genotype imputation,” Bioinformatics, vol. 31, no. 5, pp. 782–784, 2015.

S. A. Bacanu, B. Devlin, and K. Roeder, “The power of genomic control,” American Journal of Human Genetics, vol. 66, no. 6, pp. 1933–1944, 2000.

C. J. Willer, Y. Li, and G. R. Abecasis, “METAL: fast and efficient meta-analysis of genomewide association scans,” Bioinformatics, vol. 26, no. 17, pp. 2190-2191, 2010.

C. Chelala, A. Khan, and N. R. Lemoine, “SNPnexus: a web database for functional annotation of newly discovered and public domain single nucleotide polymorphisms,” Bioinformatics, vol. 25, no. 5, pp. 655–661, 2009.

L. D. Ward and M. Kellis, “HaploReg: a resource for exploring chromatin states, conservation, and regulatory motif alterations within sets of genetically linked variants,” Nucleic Acids Research, vol. 40, no. D1, pp. D930–D934, 2011.

Q. Zheng and X. J. Wang, “GOEAST: a web-based software toolkit for Gene Ontology enrichment analysis,” Nucleic Acids Research, vol. 36, Supplement 2, pp. W358–W363, 2008.

R. E. Bell, T. Golan, D. Sheinboim et al., “Enhancer methylation dynamics contribute to cancer plasticity and patient mortality,” Genome Research, vol. 26, no. 5, pp. 601–611, 2016.

S. B. Baylin and P. A. Jones, “A decade of exploring the cancer epigenome - biological and translational implications,” Nature Reviews Cancer, vol. 11, no. 10, pp. 726–734, 2011.

L. M. Boyden, J. Mao, J. Belsky et al., “High bone density due to a mutation in LDL-receptor-related protein 5,” The New England Journal of Medicine, vol. 346, no. 20, pp. 1513–1521, 2002.

M. V. Semenov and X. He, “LRP5 mutations linked to high bone mass diseases cause reduced LRP5 binding and inhibition by SOST,” The Journal of Biological Chemistry, vol. 281, no. 50, pp. 38276–38284, 2006.

J. A. Riancho, J. M. Olmos, B. Pineda et al., “Wnt receptors, bone mass, and fractures: gene-wide association analysis of LRP5 and LRP6 polymorphisms with replication,” European Journal of Endocrinology, vol. 164, no. 1, pp. 123–131, 2010.

J. B. Richards, H. F. Zheng, and T. D. Spector, “Genetics of osteoporosis from genome-wide association studies: advances and challenges,” Nature Reviews Genetics, vol. 13, no. 8, pp. 576–588, 2012.

B. Tuysuz, A. Bursali, Z. Alp, N. Suyugul, C. M. Laine, and O. Makitie, “Osteoporosis-pseudoglioma syndrome: three novel mutations in the LRP5 gene and response to bisphosphonate treatment,” Hormone Research in Pædiatrics, vol. 77, no. 2, pp. 115–120, 2012.

M. Borrell-Pages, J. C. Romero, and L. Badimon, “LRP5 negatively regulates differentiation of monocytes through abrogation of Wnt signalling,” Journal of Cellular and Molecular Medicine, vol. 18, no. 2, pp. 314–325, 2014.

M. Kimura and Y. Okano, “Molecular cloning and characterization of the human NIMA-related protein kinase 3 gene (NEK3),” Cytogenetic and Genome Research, vol. 95, no. 3-4, pp. 177–182, 2001.

S. L. Miller, G. Antico, P. N. Raghunath, J. E. Tomaszewski, and C. V. Clevenger, “Nek 3 kinase regulates prolactin-mediated cytoskeletal reorganization and motility of breast cancer cells,” Oncogene, vol. 26, no. 32, pp. 4668–4678, 2007.

J. Chang, R. H. Baloh, and J. Milbrandt, “The NIMA-family kinase Nek 3 regulates microtubule acetylation in neurons,” Journal of Cell Science, vol. 122, no. 13, pp. 2274–2282, 2009.

T. Zeller, P. Wild, S. Szymczak et al., “Genetics and beyond—the transcriptome of human monocytes and disease susceptibility,” PLoS One, vol. 5, no. 5, article e10693, 2010.

H. Takayanagi, S. Kim, T. Koga et al., “Induction and activation of the transcription factor NFATc1 (NFAT2) integrate RANKL signaling in terminal differentiation of osteoclasts,” Developmental Cell, vol. 3, no. 6, pp. 889–901, 2002.

M. Matsumoto, M. Kogawa, S. Wada et al., “Essential role of p 38 mitogen-activated protein kinase in cathepsin K gene expression during osteoclastogenesis through association of NFATc1 and PU.1,” The Journal of Biological Chemistry, vol. 279, no. 44, pp. 45969–45979, 2004.

K. Kim, J. H. Kim, J. Lee et al., “Nuclear factor of activated T cells c1 induces osteoclast-associated receptor gene expression during tumor necrosis factor-related activation-induced cytokine-mediated osteoclastogenesis,” The Journal of Biological Chemistry, vol. 280, no. 42, pp. 35209–35216, 2005.

Y. Kim, K. Sato, M. Asagiri, I. Morita, K. Soma, and H. Takayanagi, “Contribution of nuclear factor of activated T cells c1 to the transcriptional control of immunoreceptor osteoclast-associated receptor but not triggering receptor expressed by myeloid cells-2 during osteoclastogenesis,” The Journal of Biological Chemistry, vol. 280, no. 38, pp. 32905–32913, 2005.

M. M. Winslow, M. Pan, M. Starbuck et al., “Calcineurin/NFAT signaling in osteoblasts regulates bone mass,” Developmental Cell, vol. 10, no. 6, pp. 771–782, 2006.

J. H. Kim and N. Kim, “Regulation of NFATc1 in osteoclast differentiation,” Journal of Bone Metabolism, vol. 21, no. 4, pp. 233–241, 2014.

R. J. Pruim, R. P. Welch, S. Sanna et al., “LocusZoom: regional visualization of genome-wide association scan results,” Bioinformatics, vol. 26, no. 18, pp. 2336-2337, 2010.