ORIGINAL RESEARCH article

Front. Cell Dev. Biol., 31 March 2020

Sec. Epigenetics and Genome Architecture

Volume 8 - 2020 | https://doi.org/10.3389/fcell.2020.00184

Characterization of Dysregulated lncRNA-Associated ceRNA Network Reveals Novel lncRNAs With ceRNA Activity as Epigenetic Diagnostic Biomarkers for Osteoporosis Risk

  • Department of Geriatrics, The Second Affiliated Hospital of Harbin Medical University, Harbin Medical University, Harbin, China

Abstract

The altered expression of long non-coding RNAs (lncRNAs) has been implicated in the development and human diseases. However, functional roles and regulatory mechanisms of lncRNA as competing endogenous RNAs (ceRNAs) in osteoporosis and their potential clinical implication for osteoporosis risk are largely unexplored. In this study, we performed integrated analysis for paired expression profiles and regulatory relationships of dysregulated lncRNAs, mRNAs, and miRNAs based on “ceRNA hypothesis,” and constructed an osteoporosis-related dysregulated miRNA-mediated lncRNA–mRNA ceRNA network (DysCeNet) composed of 105 nodes (including eight miRNAs, 24 mRNAs, and 73 lncRNAs) and 515 edges. Functional analysis suggested that the DysCeNet was involved in known osteoporosis or bone metabolism-related biological processes and pathways. Then, we performed random forest-based feature selection for 73 lncRNAs with ceRNA activity and identified 25 of 73 lncRNAs as potential diagnostic biomarkers. A random forest-based classifier composed of 25 lncRNA biomarkers (RF-25lncRNA) was developed for predicting osteoporosis risk. Performance evaluation with the leave-one-out cross-validation (LOOCV) procedure showed that the RF-25lncRNA achieved a good performance in distinguishing high- and low-bone mineral density (BMD) subjects in different osteoporosis datasets. Our study for the first time revealed a global view of lncRNA-associated ceRNA regulation in osteoporosis and provided novel lncRNAs with ceRNA activity as candidate epigenetic diagnostic biomarkers for early detection of osteoporosis risk.

Introduction

Osteoporosis is a progressive systemic skeletal disease with low bone density and deterioration of bone architecture (; ). It is estimated that over 14 million men and women in the United States will have osteoporosis by 2020 (). People with osteoporosis substantially have an increased risk of bone fragility and fracture, leading to increase pain, disability, nursing home placement, total health care costs, and mortality (). Early detection and intervention for osteoporosis risk are an effective way to delay the development of disease and improve the quality of life of patients. Therefore, there is an urgent need to identify novel molecular biomarkers in the clinical assessment of osteoporosis risk.

With the development and application of high-throughput omics technologies, it has been shown that osteoporosis has genetic and molecular heterogeneity like many other common complex diseases, and disruption in some molecular pathways can disturb the equilibrium of bone turnover and thereby contribute to osteoporosis (). Long non-coding RNAs (lncRNAs) are a newly discovered class of non-coding RNAs (ncRNAs) that are longer than 200 bp () and have attracted much attention in recent medical studies. A large number of studies have demonstrated that lncRNAs is an important player of the genomic regulatory network and is involved in a wide variety of biological progress (; ; ; ; ). lncRNAs have been attributed to various functions in development, differentiation, and human disease by negatively or positively regulating gene expression at the transcriptional, post-transcriptional, and epigenetic levels (; ). Growing functional roles of lncRNAs has been highlighted in osteoporosis during recent studies. For example, XIST, a well-known major effector of the X-inactivation process, was recently reported to promote osteoporosis through inhibiting bone marrow mesenchymal stem cell differentiation (). found that lncRNA ZBTB40-IT1 played a critical role in bone metabolism and can be modulated by osteoporosis GWAS risk SNPs suppresses osteogenesis. It has become increasingly clear that lncRNAs can act as competing endogenous RNAs (ceRNAs) to interact other RNA molecules by competing for binding to shared microRNAs, which has been implicated in development and human disease including osteogenesis (; , ; ). However, genome-wide exploration for miRNA-mediated lncRNA-associated ceRNA mechanism in osteoporosis and their potential clinical implication for osteoporosis risk remained largely unknown.

In this study, we tried to construct a global miRNA-mediated lncRNA–mRNA ceRNA network in the development of osteoporosis by integrating paired expression profiles and regulatory relationship of lncRNAs, mRNAs, and miRNAs based on “ceRNA hypothesis,” and further to uncover novel lncRNAs with ceRNA activity as epigenetic diagnostic biomarkers for identifying people at high risk for developing osteoporosis.

Materials and Methods

Patients and Samples

Ten samples [including five high-bone mineral density (BMD) subjects and five low-BMD subjects] with corresponding transcriptome gene expression microarray data (Affymetrix Human Exon 1.0 ST Array) and epigenomic miRNA microarray data (Affymetrix Multispecies miRNA-2 Array) were obtained from the Gene Expression Omnibus (GEO) database (the accession number is GSE625891). Two other osteoporosis datasets with transcriptome gene expression microarray data were obtained from the GEO database, including GSE56814 dataset (16 high-BMD subjects and 15 low-BMD subjects) (the accession number is GSE568142) and GSE13850 dataset (10 high-BMD subjects and 10 low-BMD subjects) (the accession number is GSE138503).

Acquisition and Analysis of Expression Profiles

Raw transcriptome gene expression microarray data (CEL files) profiled on Affymetrix GeneChip Human Exon 1.0 ST Array (HuEx-1_0-st) and Affymetrix Human Genome U133A Array (HG-U133A) were obtained from the GSE63402 and GSE13850. These raw data were processed and normalized using the Robust Multichip Average (RMA) algorithm of R package “oligo” for background subtraction, quantile normalization, and summarization. Then, all probes of microarray were mapped into protein-coding genes using the R package “biomaRt” (). LncRNA expression profiles were obtained using repurposing strategy by mapping array probes into the human genome (GRCh 38) and lncRNA annotations from the GENCODE database4 () using SeqMap tool (). Human mature miRNAs were retrieved from the miRNA microarray. Finally, expression profiles of 20,068 mRNA, 7821 lncRNAs, and 1100 miRNA were obtained for further analysis.

Differential expression analysis of lncRNAs, miRNAs, and mRNAs between high- and low-BMD subjects was performed using the R package “limma” (). Those lncRNAs, miRNAs, and mRNAs with p < 0.05 were considered as differentially expressed genes. Hierarchical clustering was performed to investigate the expression patterns between high- and low-BMD subjects.

Construction of Dysregulated lncRNA-Associated ceRNA Network

The experimentally validated miRNA–mRNA and miRNA–lncRNA interaction data were collected from the TarBase database5 (). The dysregulated lncRNA-associated ceRNA network (DysCeNet) in osteoporosis was constructed based on the “ceRNA hypothesis” as follows: (i) Pearson correlation coefficient (PCC) was calculated to measure the expression correlation between differentially expressed mRNAs and lncRNAs from matched mRNA and lncRNA expression profiles. Those dysregulated lncRNA–mRNA pairs with PCC > 0.5 were considered as candidate co-dysregulated lncRNA–mRNA ceRNA crosstalk; (ii) expression correlation between differentially expressed miRNAs and differentially expressed mRNAs, and between differentially expressed miRNAs and differentially expressed lncRNAs was evaluated using PCC; (iii) for a candidate co-dysregulated lncRNA–mRNA ceRNA crosstalk, both mRNAs and lncRNAs in this lncRNA–mRNA ceRNA crosstalk are targeted and co-expressed negatively with a common miRNA; this candidate co-dysregulated lncRNA–mRNA ceRNA crosstalk was selected as dysregulated miRNA-mediated lncRNA–mRNA ceRNA crosstalk; (iv) all miRNA-mediated lncRNA–mRNA ceRNA crosstalks were integrated to form a global DysCeNet.

Identification of lncRNA Biomarkers With ceRNA Activity Using a Machine Learning Method

To identify potential lncRNA biomarkers with ceRNA activity, a random forest approach and leave-one-out cross-validation (LOOCV) were used to select optimal lncRNAs biomarkers using the R package “randomForest” and out-of-bag (OOB) error, which measure the performance of the model on the training set (; ). The OOB error will produce an unbiased estimate for the classification error, while the bagging method will decrease the chance of overfitting (). Then, a random forest-based classifier was built using the optimal lncRNA biomarkers, and a receiver operating characteristic (ROC) curve and the area under ROC curve (AUC) was used to measure the diagnostic performance of the lncRNA classifier ().

Functional Enrichment Analysis

Functional enrichment analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) for mRNAs in the DysCeNet was conducted to predict potential biological functions of lncRNAs in the DysCeNet using the R package “clusterprofiler” (). Those significantly enriched GO terms with p < 0.05 with mutually overlapping gene sets were clustered together using the Enrichment Map plugin in Cytoscape environment ().

Results

Identification of Differentially Expressed mRNAs, miRNAs, and lncRNAs Associated With Osteoporosis

To identify potential risk mRNAs, miRNAs, and lncRNAs in osteoporosis, we performed a comparative analysis for expression profiles of mRNAs, miRNAs, and lncRNAs between high- and low-BMD subjects. A total of 68 mRNAs, 11 miRNAs, and 95 lncRNAs were identified as differentially expressed (p < 0.05) in high-BMD subjects compared with low-BMD subjects (Supplementary Table S1). Among them, there were 73 unregulated genes (including six mRNAs, 10 miRNAs, and 57 lncRNAs) and 101 downregulated genes (including 62 mRNAs, one miRNA, and 38 lncRNAs) in high-BMD subjects compared with low-BMD subjects (Figure 1A). Hierarchical clustering analysis showed that expression patterns of these differentially expressed mRNAs, miRNAs, and lncRNAs were capable of distinguishing high-BMD subjects from low-BMD subjects (Figure 1B).

FIGURE 1

Construction and Analysis of Dysregulated lncRNA-Associated ceRNA Network

To construct an osteoporosis-related dysregulated miRNA-mediated lncRNA–mRNA ceRNA network, we performed an integrated analysis for paired expression profiles and regulatory relationships of dysregulated lncRNAs, mRNAs, and miRNAs of 10 samples in the GSE62589 as described in Section “Materials and Methods.” Finally, an osteoporosis-related dysregulated miRNA-mediated lncRNA–mRNA ceRNA network was constructed and was composed of 105 nodes and 515 edges (including eight miRNAs, 24 mRNAs, and 73 lncRNAs) (Figure 2A) (Supplementary Table S2). To further explore the functional implication of the DysCeNet, we performed functional enrichment analysis for mRNAs in the DysCeNet and found that mRNAs in the DysCeNet were significantly enriched in blood vessel development, NIK/NF-kappaB signaling, bone mineralization involved in bone maturation, osteoclast differentiation, and cell aging.

FIGURE 2

Identification of Potential lncRNA Biomarkers for the Osteoporosis Risk

To identify potential lncRNA biomarkers for the osteoporosis risk, we performed feature selection for lncRNAs in the DysCeNet using a random forest model. Finally, 25 lncRNAs of 73 lncRNAs in the DysCeNet were identified as potential biomarkers for the osteoporosis risk according to their discriminative power using OOB error (Table 1). Of them, 10 lncRNA biomarkers are up-regulated and 15 lncRNA biomarkers are down-regulated in high-BMD subjects compared with low-BMD subjects (Figure 3A). To test whether these 25 lncRNA biomarkers could efficiently distinguish high- and low-BMD subjects, we performed hierarchical clustering for 10 samples in the GSE62589 according to the expression pattern of 25 lncRNA biomarkers. As shown in Figure 3B, all 10 samples were classified into two clusters according to the expression pattern of seven lncRNA biomarkers with 100% accuracy. The results of hierarchical clustering demonstrated the potential of 25 lncRNAs as diagnostic biomarkers for osteoporosis risk. Therefore, a random forest-based classifier composed of 25 lncRNA biomarkers was developed. The performance of the RF-25lncRNA for predicting osteoporosis risk was evaluated in the GSE62589 dataset using the LOOCV procedure, in which nine samples were used as the training set and the remaining one was served as the test sample. Results of performance evaluation showed that the RF-25lncRNA achieves a perfect predictive performance in distinguishing high- and low-BMD subjects with an AUC of 1.0 (Figure 3C).

TABLE 1

Ensembl idGene nameLocationp-valuea
ENSG00000228290RP11-132M7.3Chr 6: 85,399,143-85,419,252 (+)0.0008
ENSG00000258610RP11-488C13.7Chr 14: 77,245,295-77,248,592 (+)0.0004
ENSG00000242009RP11-359D24.1Chr 3: 81,889,698-81,972,012 (+)0.0078
ENSG00000254562RP11-64I17.1Chr 11: 38,668,026-38,695,791 (+)0.0011
ENSG00000224239AC090044.2Chr 3: 898,807-899,774 (+)0.0063
ENSG00000279785AC008267.8Chr 7: 66,474,556-66,475,749 (+)0.0007
ENSG00000264108RP11-357H3.1Chr 18: 73,720,073-73,727,816 (−)0.0027
ENSG00000227885RP11-79N23.1Chr 6: 53,605,547-53,611,820 (−)0.0088
ENSG00000272599RP11-152N13.16Chr 10: 74,884,331-74,885,163 (−)0.0053
ENSG00000236145AC079988.3Chr 2: 121,935,903-121,940,156 (+)0.0061
ENSG00000255471RP11-736K20.5Chr 11: 86,603,256-86,636,079 (−)0.0296
ENSG00000203496AL133216.1Chr 10: 38,453,181-38,466,176 (+)0.0444
ENSG00000250893AC098869.2Chr 4: 40,426,119-40,427,585 (+)0.0299
ENSG00000278990AL132796.3Chr 14: 98,916,902-98,917,801 (−)0.0098
ENSG00000235890TSPEAR-AS1Chr 21: 44,506,807-44,516,575 (+)0.0414
ENSG00000235823OLMALINCChr 10: 100,373,099-100,454,043 (+)0.0207
ENSG00000253210AC040970.1Chr 8: 141,124,717-141,130,064 (−)0.0281
ENSG00000228705LINC00659Chr 20: 62,774,121-62,776,996 (−)0.0193
ENSG00000236054Z82246.1Chr 22: 32,583,300-32,584,204 (+)0.0097
ENSG00000251329LINC02353Chr 4: 32,351,038-32,353,220 (+)0.0171
ENSG00000234584AC019186.1Chr 2: 157,725,708-157,736,005 (−)0.0104
ENSG00000237017AC245052.4Chr 19: 54,119,511-54,125,343 (−)0.0191
ENSG00000230753ZNF341-AS1Chr 20: 33,787,373-33,811,109 (−)0.0077
ENSG00000240591AL096701.3Chr 22: 31,452,881-31,464,304 (+)0.0179
ENSG00000226649AC019118.1Chr 2: 3,156,756-3,157,797 (+)0.0087

Genomic information of 25 lncRNA biomarkers with ceRNA activity for osteoporosis risk.

ap-value was obtained from differential expression analysis.

FIGURE 3

Further Validation of lncRNA Biomarkers for the Osteoporosis Risk in Two Other Independent Datasets

To test the robustness of the lncRNA biomarkers for the osteoporosis risk, 25 lncRNA biomarkers were applied to the independent GSE56814 dataset. The RF-25lncRNA correctly classified 11 of 15 low-BMD subjects and nine of 16 high-BMD subjects, achieving an AUC of 0.733 (95% CI: 0.553–0.913) with an accuracy of 64.5% (Figure 4A).

FIGURE 4

The performance of the lncRNA biomarkers was further tested in another independent GSE13850 dataset. However, transcriptome gene expression data of the GSE13850 was profiled on the HG-U133A platform, and only two lncRNAs (RP11-488C13.7 and RP11-152N13.16) of 25 lncRNA biomarkers in the RF-25lncRNA were covered on the HG-U133A platform. When two lncRNA biomarkers (RP11-488C13.7 and RP11-152N13.16) were applied to the 20 samples of the GSE13850, two lncRNA biomarkers (RP11-488C13.7 and RP11-152N13.16) correctly classified seven of 10 low-BMD subjects and seven of 10 high-BMD subjects, achieving an AUC of 0.66 (95% CI: 0.396–0.924) with an accuracy of 70% (Figure 4B).

Discussion

Early detection and intervention for osteoporosis are crucial to prevent fragility fractures and delay disease progression. Traditional markers for osteoporosis risk are BMD, vitamin D, alkaline phosphatase, and so on (). Increasing evidence has suggested that altered molecular profiles contributed to the osteoporosis and outcome, which provided novel insights into molecular basis of osteoporosis and also highlighted the potential of molecular factors as markers for osteoporosis diagnosis and prognosis (; ; ; ). Biomarker identification has been proven to be an effective way to recognize people at high risk for developing osteoporosis and have attracted much attention in the clinical decision-making for osteoporosis management. Previous studies have focused on mRNA profiles and miRNA profiles and identified some candidate mRNA- or miRNA-based biomarkers. For example, circulating miRNAs, hsa-miR-122-5p and hsa-miR-4516, have been found to be diagnostic biomarkers for osteoporosis risk (). Another study by also identified five freely circulating miRNAs associated with osteoporosis fractures. There is increasing evidence that lncRNAs also play important roles in the pathogenesis of osteoporosis (; ; ). In this study, we obtained lncRNA profiles of osteoporosis patients by repurposing array probes on publicly available microarray data and performed genome-wide analysis for expression patterns of lncRNAs, miRNAs, and mRNAs between high- and low-BMD subjects. A total of 68 mRNAs, 11 miRNAs, and 95 lncRNAs were found to be associated with osteoporosis, which provided a candidate resource for experimental scientist for studying the molecular mechanism of osteoporosis.

Through many candidate osteoporosis-related mRNAs, miRNAs, and lncRNAs were identified in our study, regulatory relationships and mechanisms among these different types of RNA molecules in the development of osteoporosis are still unknown. It is well known that different RNA molecules can act as ceRNAs to communicate with and co-regulate each other by competing for binding to shared miRNAs (). LncRNAs have been reported as key components of the ceRNA-mediated regulatory network, and aberrant miRNA-mediated lncRNA–mRNA ceRNA crosstalk has been implicated in many human complex diseases, including cancers. However, functional roles and regulatory mechanisms of lncRNA as ceRNAs in the development of osteoporosis and their potential implication for osteoporosis are largely unexplored. To explore the ceRNA activity of lncRNAs in the development of osteoporosis, we performed integrated analysis for paired expression profiles of dysregulated 68 mRNAs, 11 miRNAs, and 95 lncRNAs and miRNA–target regulatory relationship based on “ceRNA hypothesis,” and constructed an osteoporosis-related miRNA-mediated dysregulated lncRNA–mRNA ceRNA network (DysCeNet) composed of 105 nodes and 515 edges (including eight miRNAs, 24 mRNAs, and 73 lncRNAs). Network analysis suggested that a large proportion of deregulated lncRNAs (76.8%, 73/95) in osteoporosis function as ceRNA and communicated with 24 mRNAs by competing for eight common miRNAs (Figure 2A), which implied that extensive variation in miRNAs and lncRNAs disrupted the miRNA-mediated lncRNA–mRNA ceRNA regulatory network contributing to osteoporosis at the post-transcriptional level. Functional analysis through functional enrichment analysis for mRNAs in the DysCeNet found that mRNAs as ceRNA counterparts of lncRNA biomarkers were involved in known osteoporosis or bone metabolism-related biological progression and pathways, including blood vessel development, NIK/NF-kappaB signaling, bone mineralization involved in bone maturation, osteoclast differentiation, and cell aging (Figures 2B,C). For example, blood vessels in the bone play vital roles for the formation of new bone and promoting blood vessel growth could reverse the weakening of bones and treat osteoporosis (). NIK/NF-kappaB signaling has been shown to play an important role in the positive and negative regulation of cytokine-mediated osteoclast formation and activation (). These results suggested that dysregulated expression of lncRNA ceRNAs and the resultant perturbation in miRNA-mediated lncRNA–mRNA crosstalk in the DysCeNet are involved in osteoporosis-biological processes and contributed to the osteoporosis.

A large number of studies have indicated the superior potential of lncRNAs as diagnostic and prognostic biomarkers compared with protein-coding genes due to the fact that lncRNAs were expressed in much more cell- type-, tissue-, and disease-specific patterns that are closely more associated with their function (; ). lncRNA biomarkers have been widely investigated and identified in various cancers during the past years (, , ; ). To explore the potential application of lncRNAs with ceRNA activity as diagnostic biomarkers for osteoporosis risk, we performed random forest-based feature selection for 73 lncRNAs with ceRNA activity and identified 25 of 73 lncRNAs as potential diagnostic biomarkers. To accelerate the clinical application, we also developed a random forest-based classifier (RF-25lncRNA) composed of seven lncRNA biomarkers and test the performance of the RF-25lncRNA in different osteoporosis datasets. Performance evaluation of the LOOCV procedure showed that the RF-25lncRNA achieved a good performance in distinguishing high- and low-BMD subjects in three osteoporosis datasets. These results demonstrated that seven lncRNA with ceRNA activity may become reliable and powerful epigenetic diagnostic biomarkers for early detection of osteoporosis risk.

There are some limitations to this study. First, the performance of newly identified lncRNA biomarkers was validated only in three osteoporosis datasets because of no other publicly available osteoporosis datasets. Second, the biological functions of newly identified lncRNA biomarkers are unknown, although they were found to have ceRNA activity in our study. Therefore, more laboratory and clinical researches were needed.

Statements

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: Three osteoporosis datasets were obtained from the Gene Expression Omnibus (GEO) database (GSE62589, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62589, GSE56814, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE56814, and GSE13850, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE13850).

Author contributions

YZ conceived and designed the experiments. MZ and LC analyzed the data. MZ wrote the manuscript. All authors read and approved the final manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.00184/full#supplementary-material

TABLE S1

List of differentially expressed miRNAs, mRNAs and lncRNAs between high-BMD subjects and low-BMD subjects.

TABLE S2

An osteoporosis-related dysregulated miRNA-mediated lncRNA-mRNA ceRNA network.

References

Summary

Keywords

competing endogenous RNAs, long non-coding RNAs, osteoporosis, biomark, epigenetics

Citation

Zhang M, Cheng L and Zhang Y (2020) Characterization of Dysregulated lncRNA-Associated ceRNA Network Reveals Novel lncRNAs With ceRNA Activity as Epigenetic Diagnostic Biomarkers for Osteoporosis Risk. Front. Cell Dev. Biol. 8:184. doi: 10.3389/fcell.2020.00184

Received

03 February 2020

Accepted

04 March 2020

Published

31 March 2020

Volume

8 - 2020

Edited by

Jianzhong Su, Wenzhou Medical University, China

Reviewed by

Jie Zhang, Tongji University, China; Hao Lin, University of Electronic Science and Technology of China, China

Updates

Copyright

*Correspondence: Yina Zhang,

These authors share first authorship

This article was submitted to Epigenomics and Epigenetics, a section of the journal Frontiers in Cell and Developmental Biology

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics