Genome-wide association studies (GWASs) have identified more than 500 loci for bone mineral density (BMD), but functional variants in these loci are less known. The aim of this study was to identify RNA modification-related SNPs (RNAm-SNPs) for BMD in GWAS loci. We evaluated the association of RNAm-SNPs with quantitative heel ultrasound BMD (eBMD) in 426,824 individuals, femoral neck (FN) and lumbar spine (LS) BMD in 32,961 individuals and fracture in ~1.2 million individuals. Furthermore, we performed functional enrichment, QTL and Mendelian randomization analyses to support the functionality of the identified RNAm-SNPs. We found 300 RNAm-SNPs significantly associated with BMD, including 249 m6A-, 28 m1A-, 3 m5C-, 7 m7G- and 13 A-to-I-related SNPs. m6A-SNPs in OP susceptibility genes, such as WNT4, WLS, SPTBN1, SEM1, FUBP3, LRP5 and JAG1, were identified and functional enrichment for m6A-SNPs in the eBMD GWAS dataset was detected. eQTL signals were found for nearly half of the identified RNAm-SNPs, and the affected gene expression was associated with BMD and fracture. The RNAm-SNPs were also associated with the plasma levels of proteins in cytokine-cytokine receptor interaction, PI3K-Akt signaling, NF-kappa B signaling and MAPK signaling pathways. Moreover, the plasma levels of proteins (CCL19, COL1A1, CTSB, EFNA5, IL19, INSR, KDR, LIFR, MET and PLXNB2) in these pathways were found to be associated with eBMD in Mendelian randomization analysis. This study identified functional variants and potential causal genes for BMD and fracture in GWAS loci and suggested that RNA modification may play an important role in osteoporosis.
Keywords: RNA modification; bone mineral density; gene expression; genome-wide association study.