Bone-Metabolism-Related Serum microRNAs to Diagnose Osteoporosis in Middle-Aged and Elderly Women

Abstract Objective: Postmenopausal osteoporosis (PMOP), a chronic systemic metabolic disease prevalent in middle-aged and elderly women, heavily relies on bone mineral density (BMD) measurement as the diagnostic indicator. In this study, we investigated serum microRNAs (miRNAs) as a possible screening tool for PMOP. Methods: This investigation recruited 83 eligible participants from 795 community-dwelling postmenopausal women between June 2020 and August 2021. The miRNA expression profiles in the serum of PMOP patients were evaluated via miRNA microarray (six PMOP patients and four postmenopausal women without osteoporosis (n-PMOP) as controls). Subsequently, results were verified in independent sample sets (47 PMOP patients and 26 n-PMOP controls) using quantitative real-time PCR. In addition, the target genes and main functions of the differentially expressed miRNAs were explored by bioinformatics analysis. Results: Four highly expressed miRNAs in the serum of patients (hsa-miR-144-5p, hsa-miR-506-3p, hsa-miR-8068, and hsa-miR-6851-3p) showed acceptable disease-independent discrimination performance (area under the curve range: 0.747–0.902) in the training set and verification set, outperforming traditional bone turnover markers. Among four key miRNAs, hsa-miR-144-5p is the only one that can simultaneously predict changes in BMD in lumbar spine 1–4, total hip, and femoral neck (β = −0.265, p = 0.022; β = −0.301, p = 0.005; and β = −0.324, p = 0.003, respectively). Bioinformatics analysis suggested that the differentially expressed miRNAs were targeted mainly to YY1, VIM, and YWHAE genes, which are extensively involved in bone metabolism processes. Conclusions: Bone-metabolism-related serum miRNAs, such as hsa-miR-144-5p, hsa-miR-506-3p, hsa-miR-8068, and hsa-miR-6851-3p, can be used as novel biomarkers for PMOP diagnosis independent of radiological findings and traditional bone turnover markers. Further study of these miRNAs and their target genes may provide new insights into the epigenetic regulatory mechanisms of the onset and progression of the disease.


Introduction
Postmenopausal osteoporosis (PMOP), which is caused by estrogen withdrawal, is the most frequent type of primary osteoporosis and threatens nearly half of middle-aged and elderly women worldwide [1,2]. Due to the lack of preliminary symptoms and typical features, delayed diagnosis is common in clinical practice, especially in surgical systems [3]. Fragility fracture is one of the most critical complications of PMOP, leading to high disability and mortality. In China alone, the projected cost of osteoporotic fractures may reach USD 25.4 billion by 2050 [4]. Therefore, early detection is essential for alleviating the harm of PMOP.
Bone mineral density (BMD) assessed by dual-energy X-ray absorptiometry (DXA) is widely accepted as an indispensable index for defining PMOP. However, due to differences in development between developed and developing regions in medical service supply, universal access to BMD measurement seems unlikely in the short term [5]. Meanwhile, as a type of assessment method, positive imaging features always lag behind the continuous abnormality of bone metabolism, which weakens the application value of these classical examination methods in diagnosis of the early phase of bone disease [6,7]. As a direct reflection of the changes in bone homeostasis, the re-review of bone turnover markers (BTMs) in recent years seems to provide new hope for the auxiliary diagnosis of PMOP [8,9]. However, recent studies, including our own, have confirmed that there is a limited correlation between the level of serum BTMs and BMD changes [10,11]. Regretfully, few convenient and accurate biomarkers are currently available for diagnosis in the clinic.
MicroRNA (miRNA) is a type of small noncoding single-stranded RNA with 18 to 24 nucleotides [12]. As one of the epigenetic mechanisms regulating gene expression, miR-NAs mediate the posttranscriptional gene-silencing of their target genes [13]. The potential value of miRNAs as novel biomarkers for early diagnosis, treatment, and prognosis monitoring has been well verified in diseases such as cancer, obesity, and diabetes [12,14]. Although several studies have found aberrant miRNA expression in osteoporosis-induced cells and animal models [15], owing to the complexity of the pathogenesis of PMOP in humans and the imperfect public gene expression databases, the mechanisms underlying disease occurrence and progression remain to be fully elucidated. Further study on the abundance difference of miRNAs in circulating serum under pathological conditions may provide a new method to reflect the overall state of bone metabolism and the dynamic process of BMD change. Meanwhile, changes in certain miRNAs may even carry specific information about the source tissue [16] to realize the precise diagnosis and treatment of PMOP.
The primary aim of the study was to screen differentially expressed miRNAs (DEmiR-NAs) in the serum of PMOP patients and postmenopausal-without-osteoporosis (n-PMOP) controls and to validate the feasibility of using key miRNAs as biomarkers for the clinical diagnosis of disease. As a secondary aim, this study explored the expression characteristics of key miRNAs in populations with different BMDs and at different body sites. In addition, the key target gene functions and signaling pathways related to DEmiRNAs that may be involved in PMOP onset and progression were also annotated. This study highlights a novel approach in PMOP diagnosis and provides new insights into the epigenetic regulatory mechanisms of the disease.

Participants
This study surveyed 795 community-dwelling, middle-aged and elderly female participants who were recruited from June 2020 to August 2021 at the First Affiliated Hospital of Sun Yat-sen University, Guangzhou, China. The inclusion criteria were as follows: (1) age ≥ 50 years; (2) menopausal duration ≥ 1 year; and (3) signed an informed consent form before study entry. The exclusion criteria were as follows: (1) any comorbidity that could significantly affect bone metabolism, e.g., thyroid disease, diabetes, cancer, kidney disease, or ankylosing spondylitis; (2) previous treatment with anti-osteoporosis drugs or hormones (vitamin D or/and calcium supplements were allowed), e.g., estrogen or glucocorticoids; and (3) a history of tobacco smoking or alcohol dependence within the last year. Finally, a total of 83 unrelated ethnic Han Chinese women were eligible and included in the analysis.

Anthropometric Measurements
Participants wore lightweight clothing and removed their shoes before anthropometric assessments. Height and weight were measured by the corrected mechanical weight and height scale (RGZ-120, Suhong Medical Instruments Co., Changzhou, China), with an accuracy of 0.1 cm in height and 0.1 kg in weight. The average value from 3 measurements was taken for final evaluation.
Areal BMD was measured via DXA (Lunar iDXA, GE Healthcare, Chicago, IL, USA) of the lumbar spine (LS) 1-4, total hip (TH), and femoral neck (FN). All evaluations were performed by experienced diagnostic imaging physicians. The device was calibrated daily against a standard calibration phantom according to the manufacturer's instructions. Based on prior measurements, the coefficient of variation (CV) for adult measurements is 0.8% for the LS, 0.8% for the FN, and 1.4% for the TH [10].

Biochemical and Immunological Analysis
Blood samples were collected via venipuncture, with participants having fasted overnight for at least 8 h. Whole blood was left to stand at room temperature for 30 min, and serum was then collected following centrifugation at 1200× g for 10 min at 4 • C.
The analyzers were calibrated daily before the analysis of all serum samples according to the manufacturer's protocol. The routine clinical chemistry panel, including UA, ALP, calcium, and phosphorus, was detected using an AU5800 automatic biochemistry analyzer and its corresponding reagents (Beckman Coulter, Brea, CA, USA), with intra-and interassay CVs ranging from 0.5% to 4.9%. The special clinical immunology panel, including 25(OH)D, N-MID, P1NP and β-CTX, was measured using a Cobas 6000 analyzer series and its corresponding reagents (Roche, Basel, Switzerland, CH), with intra-and interassay CVs ranging from 0.6% to 4.3%.

RNA Extraction
Total RNA was extracted from 250 µL serum using TRlzol LS reagent (Invitrogen, Life Technologies, Carlsbad, CA, USA) according to the manufacturer's protocol. The RNA quantity and quality were determined using an ND-1000 Spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA).

Microarray
miRNA microarray analysis was performed by a commercial service (Kangcheng Biotech Co., Shanghai, China). Briefly, miRNA expression profiling was performed using an Agilent Human miRNA Microarray system, 8 × 60 K array (Agilent Technologies, Santa Clara, CA, USA), containing probes for 2549 human miRNAs based on the miRBase database (http://www.mirbase.org, accessed on 1 December 2020, version 21.0). RNA labeling and hybridization on the Agilent miRNA microarray chips were performed with an Agilent Quick Amp Labeling Kit (Agilent part number [p/n]: 5190-0442) and Agilent Gene Expression Hybridization Kit (Agilent p/n: 5188-5242). The hybridization images were captured with an Agilent Microarray Scanner (Agilent p/n: G2565BA) and digitized using Agilent Feature Extraction (version 11.0.1.1).
The microarray data in this study have been deposited in the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/, accessed on 1 May 2022; Accession number: GSE201543).

qRT-PCR
To confirm the findings obtained by analyzing the miRNA profiles, qRT-PCR analysis was performed using a QuantStudio5 Real-time PCR System (Applied Biosystems, Waltham, CA, USA). cDNA was obtained from 150 ng of total RNA using M-MuLV Reverse Transcriptase (Enzymatics p/n, P7040L). The PCR amplification procedures were performed according to a previous description and repeated in triplicate [18]. The relative expression level of miRNAs was normalized to that of the internal control hsa-miR-425-5p using the 2 −∆∆Ct cycle threshold method [19]. The primer sequences for the qRT-PCR assays are listed in Supplementary Table S1.

Statistical Analysis
Statistical analyses were performed using IBM SPSS Statistics (IBM, Armonk, NY, USA; version 22.0). Data are presented as the mean ± standard deviation. Independent sample t-tests and one-way analysis of variance (ANOVA) were performed to analyze the data. Binary logistic regression analyses were performed to calculate the predicted probability of different diagnostic combinations. Receiver operating characteristic (ROC) curves were plotted to evaluate the diagnostic effects of the models. The area under the ROC curve (AUC) for different prediction models was compared using the method described by DeLong et al. [20]. To investigate associations between key miRNAs and the BMD at different body sites, Pearson's correlation and partial correlation analyses were used. Multiple linear regression models were conducted to examine the factors that influenced changes in BMD. The variance inflation factor (VIF) was calculated to assess the collinearity of independent variables, and an independent variable with VIF > 10 was considered highly collinear. Differences were considered statistically significant at p < 0.05.

Clinical Characteristics of the Participants
The characteristics of the participants in the discovery set, training set, and validation set are shown in Table 1. In the discovery set and training set, no significant differences in participant characteristics, including age, BMI, age at menopause, or menopausal duration, were observed between the PMOP patients and n-PMOP controls (p > 0.05). In the validation set, age and menopausal duration were significantly higher in the 23 PMOP patients than in the 12 n-PMOP controls (67.5 ± 8.8 vs. 58.3 ± 8.1, p = 0.005; 16.0 ± 8.5 vs. 6.3 ± 6.1, p = 0.001, respectively). The differences in BTMs and biochemical indices did not reach statistical significance between PMOP patients and n-PMOP controls in the discovery set, training set, or validation set (p > 0.05), except for calcium in the discovery set (2.10 ± 0.13 vs. 2.35 ± 0.06, p = 0.006).

Screening of Key miRNAs
The expression levels of 2549 miRNAs were measured in the discovery set. Under the criteria of adj. p < 0.05 and |Log2 Fold Change (FC)| > 1, a total of 198 DEmiRNAs were screened (Supplementary Table S2). Compared with the n-PMOP controls, 148 miR-NAs were significantly upregulated and 50 miRNAs were significantly downregulated in the PMOP patients. A volcano plot was constructed to demonstrate the profiles of the DEmiRNAs (Figure 2A), and the 50 most upregulated and downregulated miRNAs are shown in the heatmap ( Figure 2B). Among the DEmiRNAs, ten miRNAs with the highest FC in expression levels compared with the n-PMOP controls were selected as candidates for verification by qRT-PCR. The qRT-PCR results showed that the differential expression of candidate miRNAs was generally consistent with the microarray results (Supplementary Table S3).
Using the method of Delong et al. [20], the differences in the AUC among different diagnostic models were examined in the training set and validation set. The results showed that although the combination of multiple miRNAs improved the AUC to a certain extent, there was no significant difference between these AUCs in either the training set or validation set (p > 0.05).

Relative Expression Levels of Key miRNAs in Different Clinical Stages
Seventy-three participants in the training set and validation set were divided into four groups based on specific guidelines [17], and the demographics data and serum indices are shown in Supplementary Table S6. The age and menopausal duration in the severe osteoporosis group (68.7 ± 5.5 and 16.2 ± 6.5, respectively) were significantly higher than in the normal group (59.2 ± 6.3 and 5.7 ± 4.5, respectively) and osteopenia group (61.4 ± 7.8 and 10.0 ± 7.0, respectively) (p < 0.05). However, no significant difference was observed between the severe osteoporosis group and the osteoporosis group in terms of demographics (p > 0.05). ALP levels (108.2 ± 36.3) in the serum of the severe osteoporosis group were the highest among the four groups (p < 0.05). N-MID, ALP, and phosphorus were the only three serum indices that showed significant differences among the four groups (p < 0.05).
The relative expression levels of the key miRNAs in serum are shown in Figure 4. hsa-miR-144-5p, hsa-miR-506-3p, and hsa-miR-6851-3p were highly expressed in the osteoporosis group, and the differences were significant compared with the normal group and osteopenia group (p < 0.05). In addition, hsa-miR-144-5p and hsa-miR-8068 were expressed at lower levels in the osteopenia group than in the severe osteoporosis group (p < 0.05). Differences in key miRNA serum levels among subgroups may partly reflect the unique bone metabolic patterns in various stages of disease.

Correlations between Key miRNAs and BMD
According to the results of Pearson correlation analysis, hsa-miR-144-5p was significantly negatively correlated with FN BMD (p < 0.05), but there was no significant correlation with the BMD at other body sites (p > 0.05); hsa-miR-506-3p and hsa-miR-8068 were significantly associated with LS 1-4 BMD (p < 0.05), while hsa-miR-6851-3p was significantly correlated with BMD at all three body sites (p < 0.05). Additionally, age and menopausal duration were used as covariates for partial correlation analysis. Apart from hsa-miR-8068, the other candidate key miRNAs had different degrees of negative correlations with BMD at LS 1-4, TH, and FN (p < 0.05), and the correlations were further enhanced (Table 3). To further evaluate the influencing factors of the key miRNAs on BMD at different body sites, multivariable linear regression was performed using the key miRNAs as independent variables and BMD as a dependent variable. Meanwhile, variables that were significant in univariate analyses were adopted as covariates. No multicollinearity was detected among the independent variables in the regression models. Table 4 shows the results from the linear regression models. In Model 1, only hsa-miR-6851-3p was found to be significantly negatively associated with LS 1-4 BMD (β = −0.851, p = 0.007). When controlling for age, menopausal duration, N-MID, ALP, and phosphorus as covariates in Model 2, hsa-miR-6851-3p remained a strong predictor of LS 1-4 BMD (β = −0.645, p = 0.026), while hsa-miR-144-5p was determined to be the lone predictor to have significant predictive power for BMD at different body sites (LS 1-4: β = −0.265, p = 0.022; TH: β = −0.301, p = 0.005; and FN: β = −0.324, p = 0.003, respectively). No significant correlation was found between BMD and hsa-miR-506-3p or hsa-miR-8068 (p > 0.05).

Target Genes and Pathways Correlated with DEmiRNAs
On the basis of the 198 DEmiRNAs, 1945 target genes were identified using the TargetScan, miRTarBase, and miRDB databases. The top 10 miRNAs and their target genes are shown in Supplementary Table S7. In total, 1120 GO terms and 85 KEGG pathways were enriched in target genes according to the criteria adj. p < 0.05 (Supplementary Tables  S8 and S9). The top 10 enriched terms of GO and KEGG pathways are shown in Figure 5A,B. The 1945 target genes were imported into the STRING database to construct a PPI network. After excluding the isolated nodes, the final PPI network is composed of 1052 nodes and 7477 edges, as is shown in Supplementary Figure S2. The top 20 hub genes were sorted by degree (refers to the number of gene connections within the network) in descending order (Supplementary Table S10).   Ninety-four target genes for the five candidate key miRNAs were also predicted. As shown in Figure 5C, hsa-miR-340-5p and hsa-miR-144-5p cooperatively regulate the target gene CREBRF. SYNPO2 is the only predicted target gene of hsa-miR-144-5p. The function of target genes was significantly enriched in 17 GO terms (Biological Process [BP]:Cellular Component [CC]:Molecular Function [MF] = 8:4:5) and a KEGG signaling pathway ("neurotrophin signaling pathway"). For the BP classification, the top three enrichment terms were "regulation of release of sequestered calcium ion into cytosol by sarcoplasmic reticulum", "regulation of ryanodine-sensitive calcium-release channel activity", and "regulation of heart rate". Membrane, focal adhesion, and cytosol were the most highly enriched CC terms. In the MF category, "protein binding", "Smad binding", and "poly(A) RNA binding" were significantly enriched ( Figure 5D). For visualization, the PPI network of target genes is presented in Figure 5E. YY1, VIM, and YWHAE were the top three hub genes, with higher interaction levels ( Table 5).

Gene
Degree Gene Degree

Discussion
Since the first discovery of miRNAs in Caenorhabditis elegans in 1993, these molecules have been widely confirmed to exist in more than 12 types of mammalian body fluids, including serum [12,21]. The high degree of conservation [22], detectability [23], and specific spatiotemporal expression [14] of serum miRNAs make them promising biomarkers for liquid biopsy. As important regulators of gene expression, an increasing number of studies have confirmed that miRNAs are involved in regulating the pathological progression of PMOP [24][25][26]. However, because of the absence of clinical validation, none of them has been widely recommended for clinical diagnostic purposes. In this study, four key miRNAs were screened and validated in independent populations. The results revealed that hsa-miR-144-5p, hsa-miR-506-3p, hsa-miR-8068, and hsa-miR-6851-3p are potential independent biomarkers for distinguishing PMOP patients from n-PMOP controls. This finding provides additional information for PMOP diagnosis, independent of radiological findings and BTMs.
Several small clinical studies have preliminarily explored the potential of abnormally expressed circulating miRNAs as diagnostic biomarkers of PMOP [27][28][29]. For example, miR-133a, a stimulator of osteoclastogenesis, is useful for the detection of PMOP [29]. However, the results of previous studies are inconsistent, possibly due to the different populations studied. Moreover, the diagnostic accuracy and clinical usability of certain miRNAs for disease in certain clinical settings are limited. A study by Mandourah et al. reported that hsa-miR-122-5p and hsa-miR-4516 were downregulated in blood samples and could serve as potential diagnostic markers of osteoporosis, but the diagnostic values are not yet sufficient (AUC = 0.752) [28]. In our study, four key miRNAs showed high diagnostic value, with AUCs greater than 0.7 in both the training set and validation set, and the AUCs using the combined miRNAs for diagnosis were higher than 0.9.
Currently, the functions of hsa-miR-144-5p and hsa-miR-506-3p have been reported mainly in the cancer field. High plasma levels of hsa-miR-144-5p were shown to be associated with renal cell carcinoma [30], non-small-cell lung cancer [31] or glioblastoma [32]. Recent research by Zhang et al. showed that miR-144-5p can reduce bone repair and regeneration in type 2 diabetes by suppressing the expression of Smad1 [33]. In this study, linear regression analysis revealed that hsa-miR-144-5p was independently associated with BMD changes in multiple body sites, and the association between its expression pattern and PMOP progression is promising for further longitudinal research. Overexpression of hsa-miR-506-3p was found to inhibit the proliferation, migration, and invasion of cancer cells in osteosarcoma [34], prostate cancer [35], and hepatocellular carcinoma [36]. Thus far, the function and serum expression profiles of hsa-miR-8068 and hsa-miR-6851-3p have not been extensively investigated. The results of the present study demonstrate for the first time the additional values of the above miRNAs in PMOP.
The predicted expression profiles of the target genes confirmed the close association between DEmiRNAs and bone metabolism. As reported by Jeong et al., YY1 significantly inhibited Runx2-mediated transcriptional activity of osteocalcin (OCN) and ALP promoters, and knockout of this gene enhanced the osteoblast differentiation induced by BMP2 and Runx2 [37]. At the same time, YY1 can modulate the transcriptional activity of Smad, thereby regulating cell differentiation induced by the transforming growth factor superfamily signaling pathway [38]. VIM is a type III intermediate filament protein that is expressed in mesenchymal cells [39]. Overexpression of VIM in osteoblasts inhibits osteoblast differentiation, as shown by reduced ALP activity, delayed mineralization, and reduced expression of osteoblast markers [40]. This effect may be mediated by VIM competitively binding ATF4, which is required for OCN transcription and osteoblast differentiation [41]. YWHAE belongs to the 14-3-3 protein family and is involved in the transduction of signaling pathways by binding to phosphoserine-containing proteins [42]. A study on YWHAE in exosomes released from an osteoblast/osteocyte coculture system revealed that YWHAE had a positive response to mechanical stress [43]. Rivero et al. loaded 14-3-3ε protein into a scaffold, and positive stimulation of osteogenicity was observed [44]. These reports regarding the functions of target genes predicted by DEmiRNA profiles support our finding that elevated levels of differentially expressed serum miRNAs are correlated with PMOP.
Despite its strengths as discussed above, this research had certain limitations. First, we used strict inclusion criteria to obtain relatively homogenous populations. However, a power calculation was not performed a priori, and the small sample size may have resulted in decreased statistical power. Second, despite accounting for many important confounders, we cannot exclude residual confounding by unmeasured or unknown confounders. Thus, these findings need to be validated in subsequent large sample studies. In addition, further functional validation in vitro and in vivo is required to confirm the associations between these molecules and disease.

Conclusions
This study presents new clinical evidence regarding the deregulated expression of miRNAs in the serum of PMOP patients. Our results indicate that hsa-miR-144-5p, hsa-miR-506-3p, hsa-miR-8068, and hsa-miR-6851-3p target a variety of bone-metabolism-related genes and pathways and are potential independent biomarkers for clinical diagnosis of the disease, outperforming traditional BTMs. Among them, hsa-miR-144-5p is the only key miRNA that can simultaneously predict changes in BMD in LS 1-4, TH, and FN. These findings provide not only a new method for clinicians to evaluate the changes in BMD in postmenopausal women under limited conditions but also a meaningful inspiration for in-depth study of the epigenetic regulation mechanisms underlying PMOP onset and progression.    Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The datasets supporting the conclusions of this article are included within the paper and its Supplementary Materials. All other datasets used and analyzed during the study are available from the corresponding author on reasonable request.