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 (n = 5905). By meta-analysis, 9 CpG-SNPs achieved a genome-wide significance level (p < 7.86 × 10-7) for association with lumbar spine BMD and additional 15 CpG-SNPs showed suggestive significant (p < 5.00 × 10-5) 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.