RNA Modification-Related Genetic Variants in Genomic Loci Associated with Bone Mineral Density and Fracture

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.


Introduction
Osteoporosis (OP) is a chronic disease characterized by decreased bone mass and damage to bone microstructures, which increases the risk of osteoporotic fracture. Due to the aging population, OP will become a serious public health problem threatening the health of middle-aged and elderly people. Bone mineral density (BMD) is an important indicator for diagnosing OP and predicting the risk of osteoporotic fracture. Calcium, vitamin D and antioxidants are beneficial to maintain BMD, while smoking and excessive alcohol consumption have adverse effects on BMD. Changes in physiological mechanisms, such as hormone levels, oxidative stress and cell apoptosis, also affect BMD levels.
BMD has a strong genetic component at all sites, with estimates of heritability ranging from 46% to 84% [1,2]. Genome-wide association study (GWAS) is a powerful tool to identify susceptibility genetic variants of complex diseases. In the past decade, large-scale GWASs have identified over 500 OP susceptibility loci [3][4][5][6]. However, exploring how

Determination of RNAm-SNPs for BMD
In this study, we used summary-level data from two large-scale BMD GWASs to determine the potential functionally related RNAm-SNPs [5,6]. Estrada et al. performed a meta-analysis on lumbar spine (LS) and femoral neck (FN) BMD, including 32,961 individuals of European and east Asian ancestry. They identified 56 loci associated with BMD at the genome-wide significance level (p < 5.0 × 10 −8 ) [6]. Morris et al. undertook a GWAS in 426,824 individuals from the UK Biobank and identified 518 genome-wide significance loci for BMD as estimated by quantitative ultrasound of the heel (eBMD). This study also identified 13 bone fracture loci, all of which were associated with eBMD, in~1.2 million individuals [5]. The BMD GWAS datasets are all publicly available at the GEnetic Factors for OSteoporosis Consortium (GEFOS) website http://www.gefos.org/.accessed on 21 January 2022.
To identify the RNAm-SNPs in the large number of SNPs from the BMD GWAS, we obtained a set of RNAm-SNPs from the RMVar database (http://rmvar.renlab.org/download. html accessed on 21 January 2022). The RMVar database contains 1678,126 RNAm-SNPs for the nine types of RNA modifications (m 6 A, m 6 Am, m 1 A, 2 -O-Me, m 5 C, m 5 U, m 7 G, A-to-I and pseudouridine). These RNAm SNPs are divided into low (predicted), median and high confidence levels. The RMVar database characterized 31,076 RNA modification-related disease-associated variants by integrating the RNAm-SNPs with GWAS and ClinVar data, which directly provided clues to explore the relationship between RNA modification and disease of interest. Based on the RNAm-SNP information downloaded from the RMVar database, we annotated the SNPs in the BMD GWAS. The RNAm-SNPs associated with BMD were identified (significance level was set to p < 5.0 × 10 −8 ) for further analysis.

Enrichment of RNAm-SNPs in the BMD GWAS Dataset
Among the BMD-associated SNPs, we determined if RNAm-SNPs were overrepresented compared to what would be expected by chance. We randomly sampled a set of non-RNAm-SNPs with the same number of RNAm-SNPs from the GWAS dataset for eBMD and calculated the proportion of SNPs with a p value < 5.0 × 10 −8 in this set of non-RNAm-SNPs. Such an operation was repeated 1000 times. The distribution of the proportions resulting from this procedure was used as a background, and then the proportion of RNAm-SNPs with a p value < 5.0 × 10 −8 was compared with this distribution, and a p value was reported.
In addition, we performed whole-genome analysis of RNAm SNPs by functional genome-wide association analysis (FGWAS) as a functional annotation of BMD. FGWAS is a program that uses association statistics calculated across the genome and incorporates functional genomic information into a GWAS to estimate the enrichment of GWAS hits in different annotation types [12]. Compared with existing methods for GWAS, FGWAS can substantially boost the detection power for discovering important genetic variants and the gene-environmental interactions influencing phenotypes and related functions. In addition, simulation studies show that FGWAS outperforms existing GWAS methods for searching sparse signals in an extremely large search space, while controlling for the family wise error rate.

eQTL Analysis for BMD-Associated RNAm-SNPs
RNA modification is important in gene expression regulation and mRNA stability and homeostasis, so RNAm-SNPs may be associated with mRNA levels. We are interested in whether the identified BMD-associated RNAm-SNPs affect mRNA expression levels and whether the affected gene expression was associated with BMD. We performed an eQTL analysis to identify the associations between RNAm-SNPs and gene expression levels in different types of cells and tissues. We obtained the eQTL information from a public database in the HaploReg browser (http://archive.broadinstitute.org/mammals/haploreg/ haploreg.php accessed on 13 March 2022). The HaploReg browser developed by the Broad Institute can show the effect of SNPs on expression from eQTL studies, predict pathogenic variation and predict target genes and possible mechanisms of disease-related variation through systematic data mining. We focused on the association between RNAm-SNPs and genes located in cis-acting eQTLs.
Integration of GWAS data with eQTL studies is helpful to prioritize functionally relevant genes in GWAS-identified loci. Zhu et al. proposed an MR method named "summary data-based Mendelian randomization" (SMR) with this idea. In this study, we performed an SMR analysis to evaluate the associations between gene expression levels in five related tissues (whole blood, adipose tissue, skeletal muscle, liver and ovary) and eBMD and fracture by integrating eQTL data from the GTEx project [13] with BMD GWAS data described above [5]. A set of cis-eQTL summary data across the five human tissues from the GTEx project were downloaded from http://cnsgenomics.com/software/smr/ #DataResource, accessed on 17 July 2020.
We ran SMR (version 0.712) with default parameters in a command-line program. The genotype data of HapMap r23 CEU were used as the reference panel to calculate the LD correlation for SMR analysis. The genome-wide significance threshold for the SMR analysis was set to 5.0 × 10 −6 . We further conducted the heterogeneity in dependent instruments (HEIDI) test to test the 'no horizontal pleiotropy' assumption, the basic assumption of the MR study. p HEIDI > 0.05 indicated that there was one single SNP affecting BMD and gene expression.

pQTL Analysis for the BMD-Associated RNAm-SNPs
We carried out pQTL analysis for the identified RNAm-SNPs to search for plasma proteins related to BMD. Through pQTL analysis, we identified the proteins associated with BMD (pleotropic association) to further understand the physiopathology of BMD loss. The associations between RNAm-SNPs and plasma protein levels were searched in three pQTL studies. The first pQTL study was a large-scale proteomics GWAS that quantified 539 associations between protein levels and gene variants in a German cohort and replicated over half of them in an Arab and Asian cohort [14]. The summary data are available in the pGWAS Server database (http://metabolomics.helmholtz-muenchen.de/pgwas/index. php?task=download accessed on 13 December 2018). The second pQTL study characterized the genetic architecture of the human plasma proteome in 3301 participants aged 18 years or older who were recruited from two cohorts in England. Briefly, this GWAS tested the association between 10.6 million SNPs and 2994 plasma proteins [15]. The summary data are publicly available (http://www.phpc.cam.ac.uk/ceu/proteins/, accessed on 15 January 2020). The third pQTL study was a proteomics GWAS that measured 83 plasma proteins in a cohort of 3394 subjects [16]. The summary SNP data can be downloaded from https: //zenodo.org/record/264128#.X3u5ga95uUk, accessed on 6 October 2020. In addition, we looked for associations between RNAm-SNPs and cytokines concentrations. The cytokine study studied the genetic basis for plasma levels of 41 cytokines in 8293 Finns [17]. The summary statistics of this study are available at http://www.computationalmedicine.fi/ data#Cytokine_GWAS, accessed on 3 June 2019.

MR Analysis of Proteins
Horizontal pleiotropy is pervasive in MR analysis and can distort MR tests, leading to inaccurate causal estimates, loss of statistical power, and potential false positive causal relationships. To further assess whether there were potential causal effects between proteins identified by pQTL analysis and BMD, we performed weighted median, inverse-variance weighted (IVW) MR, MR-Egger, and MR pleiotropy residual sum and outlier (MR-PRESSO) analyses. The IVW method is a weighted linear regression model that combines the ratio estimates from each IV in a meta-analysis model. The premise for using this method is that all genetic variations are valid instrumental variables [18]. To obtain a more reliable MR estimate, we conducted MR-Egger regression. The MR-Egger method was performed to account for potential pleiotropy of the genetic variants that would have influenced the outcome through pathways other than the exposure, as this could have caused bias to the analytical results [19]. The weighted median, IVW and MR-Egger analyses were performed by using the Mendelian Randomization R package. MR-PRESSO is a method that systematically detects and corrects horizontal pleiotropic outliers in MR testing through three steps: the MR-PRESSO global test, the MR-PRESSO outlier test and the MR-PRESSO distortion test. The outlying genetic variants were identified by applying this method [20]. The source code and documents for MR-PRESSO are available at https://github.com/ rondolab/MR-PRESSO, accessed on 17 July 2018. The default parameters were used for the MR-PRESSO analysis.
The requisite data (i.e., SNP rs number, β, standard error, and p value) were extracted from each of the BMD GWASs and pQTL studies mentioned above and then merged by SNP to form a plain file with seven columns (i.e., SNP rs number, β for protein, standard error for protein, p value for protein, β for BMD, standard error for BMD and p value for BMD) for the MR analysis using the R language. We sorted out the pQTLs with p values less than 5 × 10 −4 as potential instrumental variables. We harmonized the genetic association between the pQTL and BMD GWAS to ensure that they reflected the same effect allele. We then conducted LD clumping on these SNPs to obtain the independent pQTL (LD r 2 < 0.001, within 10,000 kb) for each protein. LD clumping was done using the clump_data function provided by the TwoSampleMR R package with reference to the 1000 Genomes EUR population.

Functional Enrichment Analysis
KEGG (Kyoto Encyclopedia of Genes and Genomes) is a collection of databases dealing with genomes, diseases, biological pathways, drugs and chemical materials that can help us understand the functional interpretation of genes and their products as a whole network. DAVID (https://david.ncifcrf.gov/ accessed on 13 May 2022) is an online bioinformatics tool that can provide comprehensive biological function annotation information for largescale genes and proteins. In this study, we used the DAVID online tool to gain insights into the functions of the potential causal proteins.
For m 1 A-SNPs, we identified 34 functional loss m 1 A-SNPs that were significantly associated with eBMD. These 34 m 1 A-SNPs belong to the high and medium confidence categories (Supplementary Table S1). The top signal was found for rs227584 in HROB (p = 1.70 × 10 −41 ), followed by the association between rs643892 in LRP5 and eBMD (p = 9.00 × 10 −31 ). The m 1 A-SNP rs227584 in HROB was also significantly associated with fracture (p = 6.90 × 10 −11 ). The association between rs643892 in LRP5 and LS-BMD was marginally significant (p = 8.63 × 10 −5 ).  Four functional loss m 7 G-SNPs belonging to the medium confidence category were significantly associated with eBMD (Supplementary Table S1). Six functional loss A-to-I-SNPs belonging to the high confidence category were significantly associated with eBMD (Supplementary Table S1). For m 5 C modification, three functional loss m 5 C-SNPs belonging to the high confidence category were significantly associated with eBMD, including rs2229503 in SPTBN1 (p = 1.70 × 10 −50 ), rs11247975 in SPON2 (p = 1.60 × 10 −23 ) and rs9986596 in ZKSCAN4 (p = 9.60 × 10 −9 ). The m 5 C-SNP rs2229503 in SPTBN1 is a synonymous variant and is also a m 6 A-SNP; rs11247975 in SPON2 and rs9986596 in ZKSCAN4 are missense variants.

Enrichment of RNAm-SNPs in the BMD GWAS Dataset
The proportion of m 6 A-SNPs, m 1 A-SNPs, m 7 G-SNPs, A-to-I-SNPs and m 5 C-SNPs that have GWAS p values < 5.0 × 10 −8 for eBMD was significantly greater than that of the non-RNAm-SNPs (Table 1). We could not test the remaining four types of RNA methylation because of the lack of data. With the FGWAS method, it was found that SNPs associated with (p < 5.0 × 10 −8 ) eBMD were significantly enriched with m 6 A-SNPs (log2 enrichment of 2.17, 95% CI: [1.00, 2.84]). Functional enrichments for the other types of RNA methylation in the eBMD GWAS dataset were not detected by using FGWAS software.

Gene Expression Associated with BMD
To further clarify the possible functional mechanisms underlying the identified RNAm-SNPs in association with BMD, we investigated whether they were associated with gene expression. We identified 253 BMD-associated RNAm-SNPs that showed effects in different cells or tissues, and cis-eQTL signals with corresponding local genes were found for 119 RNAm-SNPs (47.0%). Some identified RNAm-SNPs affect the expression of key OP susceptibility genes (Supplementary Table S2). Rs2229503, which is an m 6  In SMR analysis of integration of BMD GWAS with eQTL data from the GTEx project, we detected significant associations between gene expression in five related tissue types (whole blood, adipose, skeletal muscle, liver and ovary) and eBMD and fracture. A total of 142 significant pleotropic associations were detected for 23 of the genes containing RNAm-SNPs (SMR p < 5.0 × 10 −6 ) (Supplementary Table S3). The number of significant associations found in each tissue was 77 in whole blood, 37 in adipose tissue, 19 in skeletal muscle, 4 in liver and 4 in ovary. In SMR analysis, we found that the expression levels of some OP susceptibility genes were significantly associated with eBMD in the five related tissue types. The expression of FGFRL1 was associated with eBMD (p = 1.61 × 10 −8 ) and fracture (p = 1.24 × 10 −5 ) in skeletal muscle ( Figure 2); the expression levels of SPTBN1 were associated with eBMD in whole blood and adipose tissue (p = 3.15 × 10 −16 and 6.55 × 10 −10 , respectively) ( Figure 4).

Plasma Proteins Related to the RNAm-SNPs
We attempted to find plasma proteins that were related to the identified RNAm-SNPs. We found 340 significant pQTL signals (p < 5.0 × 10 −5 ) for 96 RNAm-SNPs that were significantly associated with eBMD (p < 5.0 × 10 −8 ) (Supplementary Table S4). A total of 180 proteins were detected, and these proteins were enriched in specific KEGG pathways, such as cytokine-cytokine receptor interaction (p = 6.30 × 10 −4 ), the PI3K-Akt signaling pathway (p = 8.20 × 10 −4 ), the NF-kappa B signaling pathway (p = 1.30 × 10 −3 ) and the MAPK signaling pathway (p = 2.20 × 10 −3 ) ( Figure 5). Most of the pQTL signals were in trans effect, while rs12660627, rs8898, rs6815946 and rs28379706 were associated with plasma levels of proteins encoded by their local genes (CD109, CTSB, IDUA and PLXNB2, respectively). Most of these signals were found in the INTERVAL study, and signals for proteins encoded by CTSB, IL6ST, SELP, ICAM1 and MBL2 were found in more than one pQTL study. The m 6 A-SNP rs739468 in CACFD1 was associated with plasma levels of 45 proteins. Other SNPs, such as rs41302673, rs2517719, rs200991, rs35835721, rs9986596, rs36019691 and rs56405707, were associated with plasma levels of more than 10 proteins. The top signal was the association between rs2645429 and plasma levels of CTSB, followed by the association between rs739468 and plasma levels of SELE. We also found that RNAm-SNPs were significantly associated with proteins encoded by COL1A1, DOCK9, IBSP, MSRA, PILRA and PTHLH, which are known to be important in OP.
Most of these signals were found in the INTERVAL study, and signals for proteins encoded by CTSB, IL6ST, SELP, ICAM1 and MBL2 were found in more than one pQTL study. The m 6 A-SNP rs739468 in CACFD1 was associated with plasma levels of 45 proteins. Other SNPs, such as rs41302673, rs2517719, rs200991, rs35835721, rs9986596, rs36019691 and rs56405707, were associated with plasma levels of more than 10 proteins. The top signal was the association between rs2645429 and plasma levels of CTSB, followed by the association between rs739468 and plasma levels of SELE. We also found that RNAm-SNPs were significantly associated with proteins encoded by COL1A1, DOCK9, IBSP, MSRA, PILRA and PTHLH, which are known to be important in OP.

Proteins Causally Associated with BMD
pQTL analysis showed that RNAm-SNPs were associated with plasma protein levels.
To support the functional role of RNAm-SNPs in BMD, we still need to demonstrate that the plasma proteins affected by the RNAm-SNPs were associated with BMD. We chose proteins for MR analysis from three aspects based on the findings of the pQTL analysis. First, CTSB, IL6ST, SELP, ICAM1 and MBL2 showed significant signals in more than one

Proteins Causally Associated with BMD
pQTL analysis showed that RNAm-SNPs were associated with plasma protein levels.
For the five proteins that showed significant signals in more than one pQTL study, significant associations of CTSB with eBMD were detected using data from both the KORA and INTERVAL studies (Table 2). By using data from the KORA study, the plasma level of CTSB was associated with eBMD in weighted medium (p = 3.51 × 10 −10 ), IVW (p = 5.82 × 10 −3 ) and MR-Egger (p = 4.17 × 10 −6 ) analyses. We found that the plasma level of CTSB was associated with eBMD in IVW (p = 4.41 × 10 −2 ) and MR-Egger (p = 6.09 × 10 −4 ) analyses using data from the INTERVAL study. The plasma level of MBL2 was associated with eBMD in weighted medium (p = 9.86 × 10 −12 ) and MR-Egger (p = 2.41 × 10 −3 ) analyses using data from the INTERVAL study, but the association may be due to pleiotropic effects because significant associations were not detected in MR-PRESSO analysis.
For the 31 proteins that were enriched in KEGG pathways, the plasma levels of IL19, LIFR, CCL1, CCL19, INSR, EFNA5, KDR and MET were found to be associated with eBMD in the MR analyses. For the remaining proteins, associations between plasma levels of CD109, COL1A1, DOCK9, PILRA, PLXNB2 and eBMD were found. However, the associations between plasma levels of CD109 and PILRA and eBMD may be due to pleiotropic effects because significant associations were not detected in MR-PRESSO analysis. Causal evidence was strong for the associations between plasma levels of COL1A1, IL19, KDR, MET and PLXNB2 and eBMD because the associations passed all four MR analysis tests and the intercepts of MR-Egger analyses were not significant. According to the MR analyses, the plasma levels of ten proteins (CCL19, COL1A1, CTSB, EFNA5, IL19, INSR, KDR, LIFR, MET and PLXNB2) were associated with eBMD. We performed a functional enrichment for these ten proteins and found that these proteins were enriched in KEGG pathways of the PI3K-Akt signaling pathway (COL1A1, INSR

Discussion
This study examined the association between RNAm-SNPs and BMD and showed that many BMD-associated SNPs in important genes identified by GWAS were related to RNA modification types of m 6 A, m 1 A, m 5 C, m 7 G and A-to-I. The findings indicate that RNA modification may play a role in OP, as the enrichment analysis showed that GWAS signals were significantly enriched with m 6 A-SNPs. Some of the identified RNAm-SNPs were located in well-known OP susceptibility genes and pathways. These SNPs showed cis-acting eQTL effects in relevant tissues, and some of them were found to be associated with proteins that were enriched in specific pathways. Moreover, the affected gene expression and protein levels were found to be associated with BMD. The results suggest a relationship among genetic variants, gene expression and BMD, i.e., the RNAm-SNPs affect RNA modification, which controls mRNA expression, and the altered mRNA expression or protein levels result in abnormal BMD.
Great efforts have been made to explore the relationship between genetic variations and diseases. However, distinguishing pathogenic variants from a large number of genetic variants remains a challenge. Most of the identified variants are functionally neutral, and only a few can cause disease. Increasing evidence shows that mRNA modifications play a critical role in modulating biological processes such as gene expression [21], mRNA stability [22] and homeostasis [23]. Annotating the functional effects of gene variants on RNA modification may be a valuable strategy to decipher the pathogenic mechanism of diseases. It is possible to determine the causal variants by measuring RNAm-SNPs associated with BMD. In this study, we demonstrated that searching for SNPs with specific functions in BMD-associated loci was essential for a better understanding of GWAS signals. We annotated many RNAm-SNPs related to BMD and showed that RNAm-SNPs may be involved in the pathogenesis of OP. Many of the identified BMD-associated RNAm-SNPs were involved in critical BMD genes and pathways. Well-known BMD genes, such as WNT4, WLS, SPTBN1, SEM1, FUBP3, LRP5 and JAG1, contain RNAm-SNPs. We found evidence to show that the RNAm-SNPs may have impacts on the expression of these modifiable genes at both the mRNA and protein levels, and the affected gene expression was associated with BMD. Our research demonstrated how to refine the associated signals from the RNAm-SNPs identified in the GWAS dataset. The findings indicated that RNA modification may play a role in OP, as these key BMD genes were affected by RNA modification. Thus far, it is not clear how RNA modification affects these genes and contributes to OP, and its underlying mechanisms need to be clarified.
Rs2229503 in SPTBN1, rs7536301 in WNT4, rs76324150, rs4792891 and rs17650901 in MAPT, rs7350928 and rs17574425 in KANSL1, and rs643892 and rs73516825 in LRP5 are RNAm-SNPs significantly associated with eBMD. Based on the RMVar database, rs2229503 is a m 6 A-and m 5 C-related SNP in SPTBN1. According to the HaploReg database, rs2229503 shows a cis-acting eQTL effect on the SPTBN1 gene. SPTBN1 is a candidate causal gene for BMD. Several GWASs have identified the relationship between SPTBN1 and BMD [5,6], and functional studies have also suggested its functional relevance to the regulation of bone mass, supporting SPTBN1 as a promising gene in the field of OP [24][25][26][27]. Upregulation of SPTBN1 inhibited STAT3 signaling, whereas STAT3 was reported to be negative for bone homeostasis. Zhang et al. demonstrated that STAT3-deficient mice were prone to developing OP. In this way, SPTBN1 suppressed STAT3 signaling, which led to decreased BMD and bone loss. SPTBN1 passed the SMR analysis, which indicated that the expression level of SPTBN1 was associated with BMD. Rs643892 is the m 1 A RNAm-SNP of LRP5, a component of the Wnt signaling pathway. According to HaploReg datasets, rs643892 is an eQTL SNP, showing a cis-acting eQTL effect on the LRP5 gene. There is increasing evidence regarding the key role of the LRP5 gene in regulating bone metabolism. The Wnt pathway plays a key role in bone metabolism [28]. It influences the differentiation and function of osteoblasts and osteoclasts, and its dysregulation leads to various forms of inherited bone mass disorders. RNAm-SNPs in WNT4, WNT2B and APC were associated with eBMD. RNAm-SNPs in MAPT and KANSL1 were associated with the expression levels of WNT3. In addition to these genes, the expression levels of many genes affected by the RNAm-SNPs were found to be associated with BMD in relevant tissues. These RNAm-SNPs may be noteworthy functional SNPs. Modification of these sites may affect gene expression and disturb bone metabolism.
RNAm-SNPs may also affect gene expression at the protein level associated with BMD. According to pQTL analysis, RNAm-SNPs are pQTLs in genes known to be important in OP, including COL1A1, DOCK9, IBSP, MSRA, PILRA and PTHLH. The identified proteins also pointed to biological pathways highly related to OP, e.g., cytokine-cytokine receptor interaction, the PI3K-Akt signaling pathway, the NF-kappa B signaling pathway and the MAPK signaling pathway [29,30]. Most importantly, these proteins were found to be causally associated with BMD in MR analysis. The MR analysis highlighted proteins in the PI3K-Akt signaling pathway (COL1A1, INSR, KDR, EFNA5 and MET), the Rap1 signaling pathway (INSR, KDR, EFNA5 and MET), the Ras signaling pathway (INSR, KDR, EFNA5 and MET) and the MAPK signaling pathway (INSR, KDR, EFNA5 and MET), that may be candidate targets for OP. The relationships between many of the identified proteins, such as COL1A1, INSR, KDR and MET, and bone metabolism have been widely studied. MR analysis determined the risk factors for OP, and the results suggest that genes in the PI3K-Akt signaling pathway, the Rap1 signaling pathway, the Ras signaling pathway and the MAPK signaling pathway play functional roles in OP. The findings also indicate that RNAm-SNPs and RNA modification may play roles in OP through specific pathways.
This study has some potential limitations. First, the m 6 A-SNP set was large, but data for other types of RNA modification were very rare, so few RNAm-SNPs of these types were identified. Second, the functional effects of RNAm-SNPs on BMD were not examined experimentally. Further experiments in OP-related cells are needed to test their function.

Conclusions
This study annotated many RNAm-SNPs in many BMD-associated loci in GWAS, and elucidated the relationship between the SNPs and BMD-associated gene expression and protein levels. The findings of this study increase our understanding of the associations identified in BMD GWAS. The results suggest that RNA modification may be involved in the pathogenesis of OP. The RNAm-SNPs in BMD loci were associated with gene expression, including mRNA levels and protein levels, and the gene expression was associated with BMD, indicating that these genes may be causal factors for OP. No previous study has shown the relationship between these kinds of RNA modification and BMD. Therefore, this study may add valuable clues for further understanding the functional mechanism underlying the development of hypertension.
Supplementary Materials: The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/genes13101892/s1, Figure S1: Flow chart of the study design; Table S1: RNAm-SNPs identified for eBMD and fracture; Table S2: RNAm-SNPs affected expression of key OP susceptibility genes; Table S3: Gene expression associated with BMD in SMR analysis; Table S4: RNAm-SNPs affected circulating protein levels.