The Integrative Analysis of Competitive Endogenous RNA Regulatory Networks in Osteoporosis

doi:10.21203/rs.3.rs-1345754/v1

Download PDF

Research Article

The Integrative Analysis of Competitive Endogenous RNA Regulatory Networks in Osteoporosis

https://doi.org/10.21203/rs.3.rs-1345754/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background: Osteoporosis (OP) is a common bone disease of old age resulting from the imbalance between bone resorption and bone formation. circRNAs are a class of endogenous noncoding RNAs (ncRNAs) involved in gene regulation and may play important roles in the development of OP. Here, we aimed to discover the OP‑related circRNA-miRNA-mRNA (ceRNA) network and the potential mechanisms.

Results: Six datasets were obtained from GEO databases and the differentially expressed genes, circRNAs, and miRNAs were screened out from these datasets. Then, combined with the predicted target miRNAs and target genes, a total of 7 target circRNAs, 5 target miRNAs, 38 target genes, and a ceRNA network was constructed. Then RNA-seq verification by using total RNA isolated from the femurs of normal and oophorectomy Wistar rats indicated that MFAP5, CAMK2A, and RGS4 were closely associated with osteoporosis. Function enrichment analysis indicated that the target circRNAs, miRNAs, and genes were involved in the process of MAPK cascade, hormone stimulus, cadherin binding, rRNA methyltransferase, PI3K-Akt signaling pathway, and Vitamin digestion and absorption, et.al. Then a circRNA-miRNA-hub gene subnetwork was constructed and the RT-qPCR analysis of human bone tissues from the femoral head was used to confirm that the transcription of hsa_circR_0028877, hsa_circR_0082916, DIRAS2, CAMK2A, and MAPK4 showed a significant correlation with osteogenic genes. Besides, the two axes of hsa_circR_0028877 - hsa-miR-1273f - CAMK2A and hsa_circR_0028877 - hsa-miR-1273f - DIRAS2 conformed to be closely associated with OP. Additionally, by constructing a drug-target gene network, RKI-1447, FRAX486, Hyaluronic, and Fostamatinib were identified as therapeutic options for OP.

Conclusions: Our study revealed the potential links between circRNAs, miRNAs, and mRNAs in OP, suggesting that the ceRNA mechanism may contribute to the occurrence of OP.

Osteoporosis

circRNA

ceRNA

Bioinformatics

GEO

Molecular sponge

Osteoporosis (OP) is a common geriatric, systemic and metabolic skeletal disease characterized by a general impairment of bone mass and bone strength that results in fragility fractures [1]. With the aging of the population around the world, the incidence of osteoporosis is increasing rapidly and the research on the pathogenesis of osteoporosis and related medicines has become exceedingly needed.

CircRNAs belong to a new class of endogenous noncoding RNAs (ncRNA) with a covalent closed-loop structure and are generated from the linear host genes with a unique ‘back-splicing’ process [2]. circRNAs are characterized by high tissue specificity and are always with a very low abundance in various tissues when compared to the protein-coding mRNAs [3, 4]. However, researches have shown that circRNAs are rich in miRNA binding sites (miRNA response elements, MREs), which means that circRNAs can act as sponges to combine with miRNAs and to compete with mRNAs (Competing endogenous RNA mechanisms, ceRNA) [5, 6]. Currently, several studies have shown that circRNAs play important roles in the regulation of many types of diseases via the interaction with disease-associated miRNAs such as neurodegenerative diseases [7], but the relationship between circRNA and OP was still dimness. Thus, the study of the mechanisms and roles of OP relative circRNAs will lead to new insight into basic physiology and disease progression of osteoporosis.

There have been researchers who identified that some dysregulated circRNAs could influence the development of the skeletal system. For example, Dandan Zhang et. al had shown that circRNA-vgll3 could promote the osteogenic differentiation of adipose-derived mesenchymal stem cells via modulating the miRNA-dependent integrin-α5 expression [8]; Xiqiang Xu et. al had shown that circRNA-0011269 could promote osteoporosis progression through regulation of RUNX2 by combining with miR-122 [9].

In the present study, we collected the expression profiles in osteoporotic patient samples by microarray, and the differentially expressed mRNAs (also means differentially expressed genes, DEGs), differentially expressed miRNAs (DEMs) and differentially expressed circRNAs (DECs) were identified. Then by combining the predicted target genes and target miRNAs, the ceRNA (circRNA–miRNA–mRNA) regulatory network was conducted. The process flow chart was shown in Fig. 1. To evaluate the main functions of the DECs, DEGs, and DEMs, the gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed. Besides, RNA-seq analysis were conducted to verify the differential expression of target genes. Then we established a protein-protein interaction (PPI) network of target genes and extracted hub genes from the PPI network; after these, a circRNA–miRNA-hub gene subnetwork was also constructed to better understand the pathogenesis of osteoporosis, and the RT-qPCR analysis was carried out to verify the signature of the subnetwork. Besides, we also predicted the OP-associated potential therapeutic drugs from the target genes in the ceRNA network.

In brief, this study attempts to better understand the ceRNA relative pathogenesis of OP and help facilitate the improvement of the diagnostic and therapy of OP.

2.1 Identification of differential expression of circRNAs, miRNAs, and mRNAs

The volcano plots for the DECs, DEMs, and DEGs and the heat maps for the microarray data were shown in figure 2. The basic information regarding the expression of circRNAs, miRNAs, and mRNAs in the microarray datasets were listed in table 1.

Three microarray datasets (GSE116925, GSE35956, and GSE156508) were included in our study to obtain the differently expressed mRNAs. After the analysis with GEO2R tool, a total of 988 mRNAs were found in gene chip GSE35956 (adj p<0.05 and logFC>|2|); 992 mRNAs were determined in gene chip GSE156508 (adj p<0.05 and logFC>|0.5|); and 908 mRNAs were identified in gene chip GSE116925 (adj p<0.05 and logFC>|1|). Subsequently, we integrated the filtered mRNAs of the three datasets and integrated them with a Venn diagram, and finally, the 104 mRNAs which could be found in at least two data sets were defined as the DEGs (figure 3A).

As for the miRNAs, 320 miRNAs were found in gene chip GSE74209 (adj p<0.05 and logFC>|1|) and 186 miRNAs were found in gene chip GSE93883 (adj p<0.05 and logFC>|1|); then finally the 41 overlapped miRNAs were defined as the DEMs (figure 3B).

As for the circRNAs, we got 10 differentially expressed circRNAs (DECs) from the gene chip GSE161361 (adj p<0.01 and logFC>|4|) and treated them as the target circRNAs (figure 3C). The basic information regarding the target circRNAs was listed in table 2.

2.2 Identification of the target miRNAs and target genes and construction of the circRNA-miRNA-mRNA network

By intersecting predicted target miRNAs and DEMs obtained from GEO, we confirmed 5 target miRNAs finally (figure. 4A).

The predicted mRNAs which could bind to the target miRNAs were obtained by using the miRWalk and the TargetScan website. Then by intersecting the predicted mRNAs and DEGs obtained from GEO, 38 target genes were confirmed after the duplicates were removed (figure. 4B-F).

Finally, the circRNA-miRNA-mRNA network containing circRNA-miRNA pairs and miRNA-mRNA pairs was generated by cytoscape software (figure. 4G), and the heatmap of the 38 target genes included in the network was constructed (figure. 5).

2.3 RNA-seq verification of target genes

To further analyze the relationship between target genes and osteoporosis, we established a model of osteoporosis in ovariectomized rats, and the RNA-seqwas proceeded to analyze the differentially expressed genes (figure. 6 A and B). As shown in figure 6C, the gene expression of MFAP5, CAMK2A, and RGS4 was significantly different between the control and OVX rats. Besides, GO and KEGG enrichment analysis of thedifferentially expressed genes (DEGs) indicated that the down-regulated genes in the OVX rats were primarily enriched in the PI3K-Akt signaling pathway, endocrine resistance pathway, osteoblast differentiation, and bone development, et.al (Supplement figure 2A and B), while the up-regulated genes in the OVX rats were primarily enriched in the thyroid hormone signaling pathway, and HIF-1 signaling pathway, et.al (Supplement figure 2C and D). The above GO and KEGG items were closely associated withosteoporosis.

2.4 PPI network establishment and hub genes identification of the target genes

In total, after removing nodes not linked to any other nodes, 32 nodes and 74 edges were mapped in the PPI network (figure. 7A). The Cytohubba app in Cytoscape was used to identify hub genes in the PPI network and a significant module containing 5 nodes and 10 edges was identified. These highest-scoring nodes were screened as hub genes: DIRAS2, CAMK2A, MAPK4, CDC42BPA, and RGS4 (figure. 7B). Besides, it's worth noting that the CAMK2A and RGS4 in these hub genes were also significant differences in the RNA sequencing results.

Then, a circRNA-miRNA-hub genesubnetwork based on the five hub genes was constructed (figure. 7C). The binding sites of circRNA-miRNA were shown in supplement figure 1 and table 4, and the binding sites of the miRNA-hub gene were shown in table 5.

2.5 qRT-PCR analysisof relative RNAs

To further examine the relationship between relative RNAs and osteogenesis, we analyzed the correlation between the circRNAs, miRNAs, mRNAs in the ceRNA sub-network and the osteogenesis indicator OPG, in the same 12 cases of human bone tissues. As shown in figure 8, the expression of hsa_circ_0028877, hsa_circ_0082916, hsa_circ_0030712, DIRAS2, CAMK2A, and MAPK4 was significantly correlated with osteogenesis.

Besides, with the correlation analysis of miRNA-circRNA and miRNA-mRNA, we judged that the expression of hsa_circR_0028877-hsa-miR-1273f - CAMK2A and hsa_circR_0028877-hsa-miR-1273f - DIRAS2 conformed to the common regulatory characteristics of the ceRNA mechanism.

2.6 Function enrichment analysis of target genes, target miRNAs, and DECs

The function enrichment analysis of target genes was performed by the Metascape website and ggplot2 package (version 3.3.3). As shown in figure 9A and B,the enrichment analysis of target genes was mainly enriched in the‘MAPK cascade’, ‘cellular response to hormone stimulus’, and ‘regulation of oxidoreductase activity’, et.al. which revealed that the target genes might affect the development of osteoporosis in these signal pathways.

Because the circRNAs belong to a type of novel non-coding RNA, their function was still not being well understood and the functional analysis could not be directly performed on the DECs. Thus, we hypothesized that the host genes of DECs play significant roles in the ceRNA network and the analysis of the host genes could reveal the potential mechanisms of DECs. Then the ggplot2 package was applied for further enrichment analysis. As shown in figure 9C, the enrichment analysis of host genes was mainly enriched in the items of ‘cadherin binding’, ‘rRNA methyltransferase activity’, and ‘cellular response to oxidized low-density lipoprotein’, et.al.

The enrichment analysis of the target miRNAs was performed bythe mirPath v.3 tool in the DIANA Tools website. MiRPath is a miRNA pathway analysis web-server in which the enrichment analysis of miRNAs was performed indirectly through the corresponding target genes[10]. As shown in figure 9D and E, the enrichment analysis was mainly enriched in the items of ‘PI3K-Akt signaling pathway’, ‘Vitamin digestion and absorption’, ‘Thyroid hormone signaling pathway’, and ‘cell death’, et.al.

2.7 Identification of the potential drugs and construction of drug-target gene interaction network

We uploaded 38 target genes from the ceRNA network to the Drugbank, TTD, and CLUE database and obtained 62 kinds of candidate compounds (drugs). Then the drug-target gene interaction network was constructed to visualize their interactions (figure. 10A). Besides, as shown in fig. 10B, to further screen the candidate compounds, the sub-network which contains the drugs that simultaneously target multiple target genes was constructed and 4 drugs were screened out: RKI-1447, FRAX486, Hyaluronic-acid, and Fostamatinib (figure. 10C-F).

Osteoporosis(OP) is a systemic skeletal disease, which is related to advancing age and characterized by the systemic destruction of bone mass and microarchitecturethat results in fragility fractures[11]. With the increase of the aging population, the number of osteoporosis patients and particularly the postmenopausal osteoporosis patients are increasing rapidly, which could increase the medical and socioeconomic pressure around the world[1].

CircRNAs belong to the class of noncoding RNAs and were first discovered in humans in 1986[12]. circRNAs have a complete closed‐loop structure in which the 3′ and 5′ ends are linked by a process termed “back-splicing” and are lacking a poly (A) tail. The closed‐loop structure makes circRNAs resistant to exonucleases and thus more stable than linear RNAs[13, 14]. Besides, circRNAs have the characteristics of tissue, time, and disease specificity[13, 15].

circRNAs can act as miRNA sponges through binding to the specific miRNAs, thus relieving the inhibitory effect of miRNAs on their target genes and regulating the expression of target genes. This mechanism is known as the competitive endogenous RNA (ceRNA) mechanism[16], which is a research hotspot at present.

The abnormal expression of circRNAs was investigated in many diseases, such as colorectal cancer, liver cancer, andpancreatic ductal adenocarcinoma, et.al[17-19], and was found to be involved in a variety of bone homeostases relative events, such as the osteogenesis differentiation and bone absorption[9, 20]. However, studies of osteoporosis-associated circRNAs and ceRNA mechanisms were stillrarely at present.

In this study, to investigate the potential role of circRNAs in osteoporosis, we first constructed the circRNA-miRNA-mRNA (ceRNA) regulatory networks by using bioinformatics predictions combined with OP relative differential expression data from the GEO database. This network contained the interactions of 4 circRNAs, 5 miRNAs, and 38 mRNAs. In the ceRNA network, circRNAs act as miRNA sponges and regulate the expression of miRNA target transcripts. Here, we used the circRNA microarray data GSE161361 to get the DECs (hsa_circ_0028878, hsa_circ_0030712, hsa_circ_0134120, hsa_circ_0112102, hsa_circ_0051728, hsa_circ_0075939, hsa_circ_0028877, hsa_circ_0089345, hsa_circ_0082916, and hsa_circ_0072602). Then, we used the miRNA microarray data GSE74209 and GSE93883 to get the DEMs and used the intersection of DEMs and predicted miRNAs to identify the target miRNAs. After that, the mRNA microarray data GSE35956, GSE116925, and GSE156508 were used to get the DEGs and the intersection of DEGs and predicted mRNAs obtained by miRWalk[21]and Target scan[22]were used to identify the target mRNAs.

To further elucidate the mechanism of the ceRNA network, we constructed a PPI network of target genes and screened 5 hub genes (DIRAS2, CAMK2A, MAPK4, CDC42BPA, and RGS4). Xia Yi et.al have indicated that CAMK2A involves in the regulation of adipogenesis and osteoblastogenesis from human bone marrow mesenchymal stem cells[23]. Heirani-Tabasi A et.al have indicated that MAPK4 participates in the migration of adipose-derived mesenchymal stem cells under hypoxia-mimicking agents[24]. Li L et.al have indicated that RGS4 protects against osteoporosis through the LINC_00370-miR_222_3p-RGS4 axis. However, the links of these genes with circRNAs have not yet been explored. Then, in order to pick out a few more important ceRNA networks, a sub-ceRNA network launch around the hub genes was constructed and 10 circRNA-miRNA-mRNA axes were identified (DIRAS2-miR_1273f-has_circ_0030712, CAMK2A-miR_1273f-has_circ_0030712, CDC42BPA-miR_1273f-has_circ_0030712, DIRAS2-miR_1273f-has_circ_0028877, CAMK2A-miR_1273f-has_circ_0028877, CDC42BPA-miR_1273f-has_circ_0028877, MAPK4-miR-4739-has_circ_0134120, RGS4-miR-4739-has_circ_0134120, MAPK4-miR-4739-has_circ_0082916, and RGS4-miR-4739-has_circ_0082916). Besides, RNA-seq analysisverification of the target genes indicated that the expression of MFAP5, CAMK2A, and RGS4 was significantly changed in the OVX ratsand is closely associated with osteoporosis. What’s more, the qRT-PCR analysis was also employed to verify these networks. The qRT-PCR analysisrevealed that the expression of hsa_circ_0028877, hsa_circ_0082916, hsa_circ_0030712, DIRAS2, CAMK2A, and MAPK4 was significantly correlated with osteogenesis, and the hsa_circR_0028877-hsa-miR-1273f- CAMK2A axis and hsa_circR_0028877-hsa-miR-1273f-DIRAS2 axis were two potential OP related ceRNA mechanisms.

In the enrichment analysis, the mRNAs were shown to be mainly involved in the MAPK pathway, cellular response to hormone, regulation of oxidoreductase activity, et.al. Among them, the MAPK pathway was a classical osteogenesis-related pathway[25], which could influence the differentiation of human bone marrow stromal cells into osteoblasts[26]. Many types of hormones proved closely related to osteoporosis, such as estrogen and thyroid hormones; Estrogen could reduce bone resorption in women and promote osteoblast osteogenesis through Wnt/ERaandb-catenin pathway[27], while thyroid hormone could improve osteoblast differentiation and activity[28]. As for theregulation of oxidoreductase activity, studies have shown that osteoblastferroptosis-induced osteoporosis was associated with the oxidoreductase relative pathway[29]. Consistent with the above results, our results further proofed that the above pathways and the related genes were closely related to osteoporosis. The GO and KEGG analysis of the target miRNAs showed that theywere enriched in the‘PI3K-Akt signaling pathway’, ‘Vitamin digestion and absorption’, ‘Thyroid hormone signaling pathway’ and ‘cell death’, andthese items were extensively involved inosteoporosis[28, 30-32]. On account of the lack of direct enrichment analysis methods for circRNA, the GO and KEGG analyses of the host genes were performed to further annotate the biological functions of the differentially expressed circRNAs indirectly. Interestingly, the host genes were enriched in the cadherin binding, rRNA methyltransferase activity, and cellular response to oxidized low-density lipoprotein pathway. All these pathways have been reported to be associated with osteoporosis[33-35]. These data indicate that the mRNAs, miRNAs, and circRNAs in our ceRNA networksare likely to play important roles in OP.

To make the most of the ceRNA network, we used the target genes to explore the potential compounds or drugs with latent therapeutic effects for osteoporosis and constructed a drug-protein network. By using Drugbank, TTD, and CLUE database, four potential chemicals were screened out, including RKI-1447, FRAX486, Hyaluronic acid, and Fostamatinib. RKI-1447 is a kind of Rho-kinase inhibitor with anti-invasive and antitumor activities in breast cancer[36]. FRAX486 is a kind of P21-activated kinase inhibitor and was reported to inhibit the growth of prostate stromal cells[37]. Hyaluronic acid,a main extracellular matrix component of articular cartilage, has been proved to promote the chondrogenic differentiation of human amniotic mesenchymal stem cells (hAMSCs). Fostamatinibis the first spleen tyrosine kinase (Syk) inhibitor approved for the treatment of chronic immune thrombocytopenia (ITP) in adult patients. However, no study has reported their effects on osteoporosis. Further research is needed on whether they can be used in the treatment of osteoporosis. The drug-protein network we constructed might contribute to the further study of thesechemicals.

There are several limitations of this study that should be considered. First, our study has been based on a range of online databases and bioinformatics analysis methods. Next, the specific mechanisms in the ceRNA network, including the miRNA-mRNA interaction and circRNA-miRNA interaction, were not been verified by a series of molecular biology experiments, further research needs to be performed to confirm the conclusions in this study.

In conclusion, by employing a comprehensive strategy of bioinformatics analysis and experimental verification, we constructed a circRNA-miRNA-mRNA network and found that hsa_circ_0028877, hsa_circ_0082916, hsa_circ_0030712, DIRAS2, CAMK2A, and MAPK4 were significantly correlated with osteogenesis, and the axes of hsa_circR_0028877-hsa-miR-1273f- CAMK2A and hsa_circR_0028877-hsa-miR-1273f-DIRAS2 might function as ceRNA mechanisms to exert critical roles in osteoporosis. In addition, 4 bioactive chemicals (RKI-1447, FRAX486, Hyaluronic acid, and Fostamatinib) were determined as potential therapeutic agents for osteoporosis. Although these results need to be further verified, our study provides a novel insight into the pathogenesis and therapy for osteoporosis from the circRNA-miRNA-mRNA view.

5.1 Microarray data mining and identification of DEGs, DEMs, and DECs

The microarray data used in this study were screened out from the GEO database, a functional genomic database of NCBI (http://www.ncbi.nlm.nih.gov/geo). There were three mRNA expression profiles (GSE35956, GSE116925, and GSE156508), two miRNA expression profiles (GSE74209 andGSE93883), and one circRNA expression profile (GSE161361)obtained from the GEO database and then the online tool GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to dig out the differentially expressed RNA datasets.

After the data processing of the GEO datasets by GEO2R, the selected RNAs were analyzed with the website of the Venn diagram tool (http://bioinformatics.psb.ugent.be/webtools/Venn/). Then the sum of the intersections between every two filtered data sets of the mRNA datasets was defined as the differentially expressed genes (DEGs); the intersections between the two filtered data sets of miRNA datasets were defined as the differentially expressed miRNAs (DEMs); and the analysis result of the circRNA expression profiles by GEO2R was defined as the differentially expressed circRNAs (DECs, also defined as the target circRNAs).

5.2 Identification of the target miRNAs and the target genes

By using the online tool of Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/), overlapping miRNAs between the predicted target miRNAs of DECs (Obtained in the GSE161361 dataset) and the DEMs were defined as the target miRNAs.

The miRWalk (http://mirwalk.umm.uni-heidelberg.de/) and the TargetScan (http://www.targetscan.org/vert_71/) website were usedto obtain the predicted mRNAs which could bind with the target miRNAs. Overlapping mRNAs between the predicted target genes of DEMs and the DEGs were defined as the target genes. Besides, the significance of expression differences of the target genes was verified by RNA-seq.

5.3 Construction of the ceRNA network

Based on the miRNA-mRNA and circRNA-miRNA interactions, the circRNA-miRNA-mRNA regulatory network was established and visualized by Cytoscape software. In addition, only the target miRNAs that were predicted to combine with both the DECs and the target genes were selected for the establishment of the ceRNA network. Besides, the heat map of the target genes in the ceRNA network was constructed by Heml 1.0.3.7 software (http://hemi.biocuckoo.org/down.php).

5.4 PPI network establishment and hub genes identification of the target genes

To analyze the target genes selected in the ceRNA network, a PPI network was established by the STRING (https://string-db.org) and then visualized by the Cytoscape software (version 3.7.2). Then, one plugin of the cytoscape software, Cytohubba, was employed to recognize the highly interacted hub genes. Cytohubba was used to calculate the degree, betweenness centrality, and closeness centrality of mRNAs in the PPI network. Overlapping and top‐ranking genes among the three calculation methods were defined as the hub genes. Besides, the hub gene relative ceRNA network was filtered out from the original ceRNA network.

5.5 Enrichment analysis of target genes, target miRNAs, and DECs

Metascape online tool (http://metascape.org) and the R packages ggplot2 (version 3.3.3) and clusterProfiler were used for the GO/KEGG enrichment analysis of the target genes.

The GO annotation and KEGG pathway analyses of the target miRNAs were conducted by the mirPath v.3 tool on the DIANA TOOLS website (http://snf-515788.vm.okeanos.grnet.gr/).

As for the enrichment analysis of DECs, the host genes of DECs were chosen as the research objects for the GO/KEGG enrichment analysis with the R packages ggplot2 and clusterProfiler.

5.6 Prediction of OP-associated drugs and extraction of relevant target genes in the ceRNA network

Drugbank(https://go.drugbank.com/), TTD (Therapeutic target database,http://db.idrblab.net/ttd/)[38], and CLUE database (Connectivity Map Linked User Environment,https://clue.io/)were selected to predict the target genes and the related potential drug molecules for OP treatment. The predicted results and related target genes were imported into the Cytoscape software and mapped into a drug-gene network. To further screen out more significant drug molecules, the drugs that simultaneously target multiple target genes were singled out.

5.7 Animals

Twenty female wistar rats were obtained from the Experimental Animal Center of Dalian Medical University (Dalian, China). All animal studies complied the Guidelines for the Care and Use of Laboratory Animals of the National Institutes of Health and all experimental protocols in the current study and were approved by the institutional animal care committee of Dalian Medical University (Dalian, China; Certificate of Conformity: No. SCXK 2018-0007), and all Animal experiments were performed following ARRIVE. The rats were housed in pathogen-free facilities under a 12/12 h light/dark cycle and the temperature during the experiment was maintained at around 23 ± 2°C. The rats were fed with the standard laboratory diet and drink deionized water freely. Female rats were undergoing either bilateral oophorectomy or bilateral sham operation after pentobarbital sodium anesthesia. One week later, the surviving rats were divided into the Oophorectomy group (OVX, n=10) and Control group (SHAM, n=10). After 8 weeks of normal feeding, all animals were anesthetized with pentobarbital sodium and the femurs were obtained for the next experiment.

5.8 RNA isolation and RNA sequencing

Total RNA was isolated from the femurs of 6 rats (3 with ovariectomy and 3 healthy controls) using TRIzol reagent (Sigma, St. Louis, USA), following the manufacture’s protocol. The quality and quantity of RNA were tested on a Nanodrop ND-2000 spectrophotometer (Thermo Scientific). An equivalent of 5 mg RNA was utilized as sequencing samples and the RNA-seq was performed on the Illumina HiSeq™ Xten platform by IGENE-BOOK Biotechnology Ltd. (Wuhan, www.igenebook.com).The accession of our project is AJRS2210325019.

5.9 Clinical samples and qRT-PCR analysis

The cancellous bone tissue was gathered from 12 patients who received surgical treatment and were frozen and stored in liquid nitrogen.The written informed consent of each subject was obtained for this experiment and this study was approved by the ethics institute of theFirst Affiliated Hospital of Dalian Medical University(YJ-KY-FB-2021-22) and information about the patient's age, gender, and health status was recorded. The collected cancellous bone tissue was analyzed by qRT-PCR (Quantitative Real-time PCR) analysis.

The expression of potential circRNAs, miRNAs, and genes in these clinical samples was detected by the method of qRT-PCR. Reverse transcription of cDNA was conducted usingEasyScript®One-Step gDNA Removal and cDNA Synthesis SuperMix reagent kit (Transgen, Beijing, China). Quantitative Real-Time PCR was performed byTransStart®Top Green qPCR SuperMix (+Dye II) (Transgen, Beijing, China) on a qPCR instrument (Bio-rad, USA). The relative level of the target gene was calculated using the 2^-^∆^CT method. The Primers used in this study were synthesized and purchased from RiboBio Co. Ltd (Guangzhou, China),the sequences of primers were shown in table 3.

5.10 Statistical analysis

Expression correlation for RNA-RNA interactions in the qRT-PCR analysis part was performed by Person correlation coefficient using Graphpad Prism software (version 7). Results with P-value < 0.05 were considered statistically significant.

ceRNA	Endogenous noncoding RNA
OP	Osteoporosis
adj p	Adjusted p value
CLUE	Connectivity Map Linked User Environment
DECs	Differentially expressed circRNAs
DEGs	Differentially expressed mRNAs
DEMs	Differentially expressed miRNAs
FC	Fold change
GO	Gene Ontology
KEGG	Kyoto Encyclopedia of Genes and Genomes
logFC	log fold change
MREs	miRNA response elements
ncRNA	Endogenous noncoding RNA
OVX	Oophorectomy
PPI	Protein-Protein Interaction
RNA-seq	RNA-sequencing
TTD	Therapeutic target database

Availability of data and materials

The datasets GSE74209, GSE161361, GSE116925, GSE35956, and GSE156508 for this study can be found in the Gene Expression Omnibus [https://www.ncbi.nlm.nih.gov/geo/]. The sequencing data of this study can be found in the supplementary material.

Ethics Statement

All the tissues samples were collected with written informed consent following the Declaration of Helsinki and with the approval of the First Affiliated Hospital of Dalian Medical University.

Consent for publication

Not applicable.

Funding

This work was supported by the National Natural Science Foundation of China [82073851]; and Important Science Fund of Science and Technology Bureau of Liaoning Province [2020JH2/10300056].

Authors’ Contributions

HL and JD designed the overall idea of this study, conceived the experiments, analyzed the data, prepared the figures and tables, and authored the drafts of the manuscript. CW, QM, JW, and YJ collected the data from the GEO datasets and performed the experiments. HS and ML supervised this study and reviewed the drafts of the manuscript. All authors read and approved the final draft.

Competing interests

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.

Acknowledgments

The author would like to thank the members of the research group for their hard work and the guidance of Professor Huijun Sun.

Rachner TD, Khosla S, Hofbauer LC: Osteoporosis: now and the future. The Lancet2011, 377(9773):1276-1287.
Dragomir M, Calin G: Corrigendum: Circular RNAs in Cancer - Lessons Learned From microRNAs. #N/A2018, 8:307.
Wang J, Zhu S, Meng N, He Y, Lu R, Yan GR: ncRNA-Encoded Peptides or Proteins and Cancer. Mol Ther2019, 27(10):1718-1725.
Anastasiadou E, Jacob LS, Slack FJ: Non-coding RNA networks in cancer. #N/A2018, 18(1):5-18.
Kristensen L, Andersen M, Stagsted L, Ebbesen K, Hansen T, Kjems J: The biogenesis, biology and characterization of circular RNAs. Nature reviews Genetics2019, 20(11):675-691.
Lei K, Bai H, Wei Z, Xie C, Wang J, Li J, Chen Q: The mechanism and function of circular RNAs in human diseases. Exp Cell Res2018, 368(2):147-158.
D'Ambra E, Capauto D, Morlando M: Exploring the Regulatory Role of Circular RNAs in Neurodegenerative Disorders. Int J Mol Sci2019, 20(21).
Zhang D, Ni N, Wang Y, Tang Z, Gao H, Ju Y, Sun N, He X, Gu P, Fan X: CircRNA-vgll3 promotes osteogenic differentiation of adipose-derived mesenchymal stem cells via modulating miRNA-dependent integrin alpha5 expression. Cell Death Differ2021, 28(1):283-302.
Xu X, Chen Y, Tan B, Wang D, Yuan Z, Wang F: Circular RNA circ_0011269 sponges miR-122 to regulate RUNX2 expression and promotes osteoporosis progression. J Cell Biochem2020.
Vlachos I, Kostoulas N, Vergoulis T, Georgakilas G, Reczko M, Maragkakis M, Paraskevopoulou M, Prionidis K, Dalamagas T, Hatzigeorgiou A: DIANA miRPath v.2.0: investigating the combinatorial effect of microRNAs in pathways. Nucleic Acids Res2012, 40:W498-504.
Yuan X, Bi Y, Yan Z, Pu W, Li Y, Zhou K: Psoralen and Isopsoralen Ameliorate Sex Hormone Deficiency-Induced Osteoporosis in Female and Male Mice. Biomed Res Int2016, 2016:6869452.
Kos A, Dijkema R, Arnberg AC, van der Meide PH, Schellekens H: The hepatitis delta (δ) virus possesses a circular RNA. Nature1986, 323(6088):558-560.
Petkovic S, Müller S: RNA circularization strategies in vivo and in vitro. Nucleic Acids Res2015, 43(4):2454-2465.
Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H: Circular RNA: A new star of noncoding RNAs. Cancer Lett2015, 365(2):141-148.
Rybak-Wolf A, Stottmeister C, Glažar P, Jens M, Pino N, Giusti S, Hanan M, Behm M, Bartok O, Ashwal-Fluss Ret al: Circular RNAs in the Mammalian Brain Are Highly Abundant, Conserved, and Dynamically Expressed. Mol Cell2015, 58(5):870-885.
Qi X, Zhang D, Wu N, Xiao J, Wang X, Ma W: ceRNA in cancer: possible functions and clinical implications. J Med Genet2015, 52(10):710-718.
Ding B, Yao M, Fan W, Lou W: Whole-transcriptome analysis reveals a potential hsa_circ_0001955/hsa_circ_0000977-mediated miRNA-mRNA regulatory sub-network in colorectal cancer. Aging2020, 12(6):5259-5279.
Shen H, Liu B, Xu J, Zhang B, Wang Y, Shi L, Cai X: Circular RNAs: characteristics, biogenesis, mechanisms and functions in liver cancer. J Hematol Oncol2021, 14(1):134.
Rong Z, Shi S, Tan Z, Xu J, Meng Q, Hua J, Liu J, Zhang B, Wang W, Yu Xet al: Circular RNA CircEYA3 induces energy production to promote pancreatic ductal adenocarcinoma progression through the miR-1294/c-Myc axis. Mol Cancer2021, 20(1):106.
Chen X, Ouyang Z, Shen Y, Liu B, Zhang Q, Wan L, Yin Z, Zhu W, Li S, Peng D: CircRNA_28313/miR-195a/CSF1 axis modulates osteoclast differentiation to affect OVX-induced bone absorption in mice. RNA Biol2019, 16(9):1249-1262.
Sticht C, De La Torre C, Parveen A, Gretz N: miRWalk: An online resource for prediction of microRNA binding sites. PLoS One2018, 13(10):e0206239.
Agarwal V, Bell G, Nam J, Bartel D: Predicting effective microRNA target sites in mammalian mRNAs. #N/A2015, 4.
Yi X, Wu P, Liu J, He S, Gong Y, Xiong J, Xu X, Li W: Candidate kinases for adipogenesis and osteoblastogenesis from human bone marrow mesenchymal stem cells. Mol Omics2021.
Heirani-Tabasi A, Mirahmadi M, Mishan M, Naderi-Meshkin H, Toosi S, Matin M, Bidkhori H, Bahrami A: Comparison the effects of hypoxia-mimicking agents on migration-related signaling pathways in mesenchymal stem cells. Cell Tissue Banking2020, 21(4):643-653.
Greenblatt MB, Shim JH, Zou W, Sitara D, Schweitzer M, Hu D, Lotinun S, Sano Y, Baron R, Park JMet al: The p38 MAPK pathway is essential for skeletogenesis and bone homeostasis in mice. J Clin Invest2010, 120(7):2457-2473.
Wang Z, Liu Q, Liu C, Tan W, Tang M, Zhou X, Sun T, Deng Y: Mg in β-TCP/Mg-Zn composite enhances the differentiation of human bone marrow stromal cells into osteoblasts through MAPK-regulated Runx2/Osx. J Cell Physiol2020, 235(6):5182-5191.
Kim RY, Yang HJ, Song YM, Kim IS, Hwang SJ: Estrogen Modulates Bone Morphogenetic Protein-Induced Sclerostin Expression Through the Wnt Signaling Pathway. Tissue Eng Part A2015, 21(13-14):2076-2088.
Lademann F, Tsourdi E, Rijntjes E, Kohrle J, Hofbauer LC, Heuer H, Rauner M: Lack of the Thyroid Hormone Transporter Mct8 in Osteoblast and Osteoclast Progenitors Increases Trabecular Bone in Male Mice. Thyroid2020, 30(2):329-342.
Yang RZ, Xu WN, Zheng HL, Zheng XF, Li B, Jiang LS, Jiang SD: Exosomes derived from vascular endothelial cells antagonize glucocorticoid-induced osteoporosis by inhibiting ferritinophagy with resultant limited ferroptosis of osteoblasts. J Cell Physiol2021.
Fang Y, Xue Z, Zhao L, Yang X, Yang Y, Zhou X, Feng S, Chen K: Calycosin stimulates the osteogenic differentiation of rat calvarial osteoblasts by activating the IGF1R/PI3K/Akt signaling pathway. Cell Biol Int2019, 43(3):323-332.
Ahmadieh H, Arabi A: Vitamins and bone health: beyond calcium and vitamin D. Nutr Rev2011, 69(10):584-598.
Tian Q, Qin B, Gu Y, Zhou L, Chen S, Zhang S, Zhang S, Han Q, Liu Y, Wu X: ROS-Mediated Necroptosis Is Involved in Iron Overload-Induced Osteoblastic Cell Death. Oxid Med Cell Longev2020, 2020:1295382.
Yang Y, Kabir M, Park J, Park J, Kang J, Ju S, Shin Y, Lee S, Lee J, Kim Set al: Plasma proteomic profiling of young and old mice reveals cadherin-13 prevents age-related bone loss. Aging2020, 12(9):8652-8668.
Yu J, Xiao M, Ren G: Long non-coding RNA XIST promotes osteoporosis by inhibiting the differentiation of bone marrow mesenchymal stem cell by sponging miR-29b-3p that suppresses nicotinamide N-methyltransferase. #N/A2021, 12(1):6057-6069.
Mazière C, Salle V, Gomila C, Mazière J: Oxidized low density lipoprotein enhanced RANKL expression in human osteoblast-like cells. Involvement of ERK, NFkappaB and NFAT. Biochimica et biophysica acta2013, 1832(10):1756-1764.
Patel R, Forinash K, Pireddu R, Sun Y, Sun N, Martin M, Schönbrunn E, Lawrence N, Sebti S: RKI-1447 is a potent inhibitor of the Rho-associated ROCK kinases with anti-invasive and antitumor activities in breast cancer. Cancer Res2012, 72(19):5025-5034.
Wang Y, Gratzke C, Tamalunas A, Wiemer N, Ciotkowska A, Rutz B, Waidelich R, Strittmatter F, Liu C, Stief Cet al: P21-Activated Kinase Inhibitors FRAX486 and IPA3: Inhibition of Prostate Stromal Cell Growth and Effects on Smooth Muscle Contraction in the Human Prostate. PLoS One2016, 11(4):e0153312.
Wang Y, Zhang S, Li F, Zhou Y, Zhang Y, Wang Z, Zhang R, Zhu J, Ren Y, Tan Yet al: Therapeutic target database 2020: enriched resource for facilitating research and early development of targeted therapeutics. Nucleic Acids Res2020, 48:D1031-D1041.

Tables are only available as a download in the Supplemental Files section.

No competing interests reported.

table1BasicinformationofthemicroarraydatasetsfromGEO.xlsx
Table 1 Basic information of the microarray datasets from GEO
table2Basiccharacteristicsofthe10differentiallyexpressedcircRNAs.xlsx
Table 2 Basic information of the 10 target circRNAs
table3primer.xlsx
Table 3 Primers used in the RT-qPCR
table4PredictedcircRNAmiRNAbindingactivity.xlsx
Table 4 Predicted circRNA–miRNA binding sites and activity
table5targetscan.xlsx
Table 5 Predicted hub gene–miRNA binding sites and activity
supplementfig1.tif
Supplement figure 1 Predicted binding sites between circRNA and miRNA. The points on the ring represent reverse splicing sites.
supplementfig2.tif
Supplement figure 2 Enrichment analysis of RNA-seq results. (A) (B) GO and KEGG analysis of down-regulated genes in the OVX rats. (C) (D) GO and KEGG analysis of up-regulated genes in the OVX rats.

Download PDF

Editorial decision: Major revision
14 Mar, 2022
Reviews received at journal
14 Mar, 2022
Reviews received at journal
16 Feb, 2022
Reviewers agreed at journal
15 Feb, 2022
Reviewers invited by journal
15 Feb, 2022
Editor assigned by journal
15 Feb, 2022
Editor invited by journal
13 Feb, 2022
Submission checks completed at journal
13 Feb, 2022
First submitted to journal
10 Feb, 2022

You are reading this latest preprint version

The Integrative Analysis of Competitive Endogenous RNA Regulatory Networks in Osteoporosis

Status:

Version 1

Abstract

Figures

1 Background

2 Results

2.1 Identification of differential expression of circRNAs, miRNAs, and mRNAs

2.2 Identification of the target miRNAs and target genes and construction of the circRNA-miRNA-mRNA network

2.3 RNA-seq verification of target genes

2.4 PPI network establishment and hub genes identification of the target genes

2.5 qRT-PCR analysisof relative RNAs

2.6 Function enrichment analysis of target genes, target miRNAs, and DECs

2.7 Identification of the potential drugs and construction of drug-target gene interaction network

3 Discussion

4 Conclusion

5 Methods

5.1 Microarray data mining and identification of DEGs, DEMs, and DECs

5.2 Identification of the target miRNAs and the target genes

5.3 Construction of the ceRNA network

5.4 PPI network establishment and hub genes identification of the target genes

5.5 Enrichment analysis of target genes, target miRNAs, and DECs

5.6 Prediction of OP-associated drugs and extraction of relevant target genes in the ceRNA network

5.7 Animals

5.8 RNA isolation and RNA sequencing

5.9 Clinical samples and qRT-PCR analysis

5.10 Statistical analysis

Abbreviations

Declarations

Availability of data and materials

Ethics Statement

Consent for publication

Funding

Authors’ Contributions

Competing interests

Acknowledgments

References

Tables

Additional Declarations

Supplementary Files

Status:

Version 1