Identification and experimental validation of key m6A modification regulators as potential biomarkers of osteoporosis

Abstract

Osteoporosis (OP) is a severe systemic bone metabolic disease that occurs worldwide. During the coronavirus pandemic, prioritization of urgent services and delay of elective care attenuated routine screening and monitoring of OP patients. There is an urgent need for novel and effective screening diagnostic biomarkers that require minimal technical and time investments. Several studies have indicated that N6-methyladenosine (m6A) regulators play essential roles in metabolic diseases, including OP. The aim of this study was to identify key m6A regulators as biomarkers of OP through gene expression data analysis and experimental verification. GSE56815 dataset was served as the training dataset for 40 women with high bone mineral density (BMD) and 40 women with low BMD. The expression levels of 14 major m6A regulators were analyzed to screen for differentially expressed m6A regulators in the two groups. The impact of m6A modification on bone metabolism microenvironment characteristics was explored, including osteoblast-related and osteoclast-related gene sets. Most m6A regulators and bone metabolism-related gene sets were dysregulated in the low-BMD samples, and their relationship was also tightly linked. In addition, consensus cluster analysis was performed, and two distinct m6A modification patterns were identified in the low-BMD samples. Subsequently, by univariate and multivariate logistic regression analyses, we identified four key m6A regulators, namely, METTL16, CBLL1, FTO, and YTHDF2. We built a diagnostic model based on the four m6A regulators. CBLL1 and YTHDF2 were protective factors, whereas METTL16 and FTO were risk factors, and the ROC curve and test dataset validated that this model had moderate accuracy in distinguishing high- and low-BMD samples. Furthermore, a regulatory network was constructed of the four hub m6A regulators and 26 m6A target bone metabolism-related genes, which enhanced our understanding of the regulatory mechanisms of m6A modification in OP. Finally, the expression of the four key m6A regulators was validated in vivo and in vitro, which is consistent with the bioinformatic analysis results. Our findings identified four key m6A regulators that are essential for bone metabolism and have specific diagnostic value in OP. These modules could be used as biomarkers of OP in the future.

Keywords: M6A; RNA modification; biomarker; bone metabolism; osteoclast; osteoporosis.

FIGURE 1

Flowchart and analysis strategy used in this study.

FIGURE 2

Standardization of gene expression. (A) Box plots of the gene expression data after normalization. (B) The uniform manifold approximation and projection (UMAP) plot and (C) principal component analysis (PCA) plot show the differences in gene expression between the two groups. The blue points represent the high bone mineral density (BMD) data, and the red points represent the low BMD data.

FIGURE 3

Expression and correlation of 14 m6A regulators in osteoporosis. (A) Protein–protein interactions of 14 m6A regulators. (B) Expression correlations of the 14 m6A regulators in all samples. The depth of the color block represents the level of the correlation coefficient, and * denotes the significance of the statistical analysis. The most correlated gene pair was YTHDF2 and RBM15, the expression status of which is presented in the scatter plot in the right panel. (C,D) The box plot and heatmap plot show the summary of 14 m6A regulators between the high- and low-BMD groups, and six m6A regulators (METTL3, METTL16, CBLL1, FTO, YTHDF2, and HNRNPC) were significantly dysregulated.

FIGURE 4

Relationship between bone metabolism-related gene sets and m6A regulators. (A) Differences in abundance and activity of bone metabolism-related gene sets in the high- and low-BMD groups. (B) Expression correlations of these bone metabolism-related gene sets and m6A regulators in all samples. Significantly, the RBM15-module pair was the most negatively correlated, the expression status of which is presented in the scatter plot in the upper right panel, while the RBM15B-multinuclear osteoclast differentiation pair was the most positively correlated with expression status presented in the scatter plot in the lower right panel.

FIGURE 5

Unsupervised clustering analysis based on the 14 m6A regulators. (A,B) Consensus clustering cumulative distribution function (CDF) and the relative area under the CDF curve for k = 2–10. According to the recommendations for selecting the number of clusters, the number of clusters with the highest average consistency was k = 2. (C) The heatmap shows the consensus matrix for the optimal k = 2. (D) The PCA plot confirmed the striking difference between the two m6A modification patterns.

FIGURE 6

Expression of the m6A regulators and bone metabolism characteristics between the two m6A modification patterns. (A) Expression differences of 14 m6A regulators in the two m6A modification patterns. (B) Differences in the abundance and activity of bone metabolism-related gene sets in the two m6A modification patterns.

FIGURE 7

Construction of a diagnostic model of osteoporosis. (A) Univariate logistic regression analysis was conducted to identify the critical m6A regulators, indicating that five m6A regulators were significant for osteoporosis (METTL16, CBLL1, YTHDF2, HNRNPC, and FTO). (B) Multivariate logistic regression was employed to identify the independent modules, and four vital m6A regulators were obtained for the diagnostic model (METTL16, CBLL1, YTHDF2, and FTO). (C) The risk score was calculated based on the expression of the four vital m6A regulators, and the median risk score (–0.366) was used as the cutoff point. All samples were divided into two groups: the high-risk group and the low-risk group. (D,E) The sensitivity and specificity of the diagnostic model in the training dataset (D) and test dataset (E) were determined by receiver operating characteristic (ROC) curves, and the area under the curve was calculated.

FIGURE 8

Creation of a regulatory network of m6A regulators-m6A target genes. (A) m6A target genes were obtained from M6A2Target, and 306 genes were potentially coregulated with the four m6A regulators. (B) Twenty-six m6A target genes were also closely related to bone metabolism. (C) KEGG pathway analysis showed that these 26 genes were mainly enriched in parathyroid hormone synthesis, secretion, action, human papillomavirus infection, and the P13K-AKT signaling pathway. (D) The GO analysis indicated that the biological processes of these genes were mainly enriched in ossification, regulation of ossification, connective tissue development, and osteoblast differentiation. (E) A regulatory network was built with four hub m6A regulators and 26 m6A target genes.

FIGURE 9

Validation of the expression of the key m6A regulators in vitro and in vivo. (A,B) Tartrate-resistant acid phosphatase (TRAP) staining and TRAP⁺ multinuclear cells counting of RAW 264.7 cells with or without nuclear factor (NF)-κB (RANKL) stimulation. Scale bar, 100 μm. (C) The m6A level in total RNA isolated from RAW264.7 cells during the osteoclast differentiation. (D) The expression of Mettl16, Fto, Cbll1, and Ythdf2 in RAW264.7 cells was detected by real-time PCR after cultured with RANKL for 4 days (E,F) Western blotting and quantification of METTL16, FTO, CBLL1 and YTHDF2 in RAW264.7 cells after cultured in RANKL. (G,H) TRAP staining and TRAP⁺ multinuclear cells counting of bone marrow-derived macrophages (BMMs) with or without RANKL and macrophage colony-stimulating factor (MCSF) stimulation. Scale bar, 100 μm. (I) The m6A level in total RNA isolated from BMMs during the osteoclast differentiation. (J) The expression of Mettl16, Fto, Cbll1, and Ythdf2 in BMMs was detected by real-time PCR after cultured with RANKL and MCSF for 4 days (K,L) Western blotting and quantification of METTL16, FTO, CBLL1, and YTHDF2 in BMMs after cultured in RANKL and MCSF. (M) A schematic diagram shows how peripheral blood monocytes (PBMCs) were obtained from the ovariectomized (OVX) and sham mice. (N) Representative images of Hematoxylin and eosin (HE) staining of mouse femurs showing the reduction of bone formation in the OVX mice relative to the sham-control counterparts. (O) The m6A level in total RNA isolated from PBMCs of the OVX and sham mice. (P) The mRNA expression level of Mettl16, Fto, Cbll1, and Ythdf2 in PBMCs of the OVX and sham groups. Compared with the sham group.