Bioinformatics analysis of microRNA linked to ubiquitin proteasome system in traumatic osteonecrosis of the femoral head.

MicroRNAs (miRNAs) have been suggested to act critical roles in the pathophysiology of traumatic osteonecrosis of the femoral head (TONFH). Unfortunately, their roles in the development of TONFH are still ambiguous. The purpose of this study is to identify promising miRNA biomarkers in traumatic osteonecrosis development.We conducted a comprehensive bioinformatics analysis using microarray datasets downloaded from the Gene Expression Omnibus database, and compared the expression of miRNAs in the serum of TONFH patients with controls. Next, we performed target prediction, function enrichment analysis, and protein-protein interaction network analysis based on differentially expressed (DE) miRNAs.We identified 26 DE miRNAs that may contribute to the pathophysiology of TONFH. The miRNAs were linked to ubiquitin proteasome system including conjugating protein ligase activity, ubiquitin-protein ligase activity and ubiquitin mediated proteolysis 5 pathway, and we exposed miR-181a-5p and miR-140-5p as promising biomarkers in TONFH.A predicting model consisting of 5 miRNAs may help discriminating high-risk patients who might develop TONFH after femur neck fracture. Among DE miRNAs, MiR-181a-5p and miR-140-5p may contribute to the development femoral head osteonecrosis after femur neck fracture via ubiquitin proteasome system.

Introduction

Osteonecrosis of the femoral head (ONFH), also named vascular or ischemic necrosis of the femoral head, is a multifactorial condition with substantial morbidity and mortality. It subsequently progresses to femoral head collapse without effective treatment, thus necessitating total hip arthroplasty. It has been estimated that 20,000 to 30,000 new cases are affected with this prevalent and disabling disorder annually in the United States.[ Moreover, around 10 percent of the 500,000 total hip arthroplasty conducted each year in the United States involve cases with ONFH.[ ONFH is generally categorized as traumatic ONFH (TONFH) and non-traumatic ONFH (NONFH). Of these, TONFH is seen commonly in patients following dislocations of the hip, femoral neck fractures, or other trauma causes. Although growing evidence has alluded to several pathogenetic mechanisms for osteonecrosis development, the exact pathogenetic mechanism has yet to be pinpointed.

MicroRNAs (miRNAs) have recently been indicated to act important roles in the pathophysiology of TONFH. MiRNAs are a set of endogenous, small, non-coding RNA molecules with an average of 23 nucleotides in length.[ They are responsible for negative regulation of gene expression after transcription by suppressing the expression of target messenger RNAs (mRNAs).[ So far, a good mount of evidence has indicated that the dysfunction of miRNAs was related to the regulation of almost all biological processes like cell proliferation, differentiation, apoptosis, cell cycle, and DNA repair,[ and also in numerous types of diseases, such as cancers, cardiovascular diseases, autoimmune diseases, as well as orthopedics diseases.[ Furthermore, accumulating evidence has suggested that miRNAs act paramount roles in the pathophysiology of ONFH. Yuan et al[ identified a total of 12 up-regulated miRNAs (eg, miR-181c-3p, miR-146a-5p, miR-99a-3p, miR-181a-3p, miR-3064-5p) and 5 down-regulated (eg, miR-132-3p, miR-212-5p, miR-629-3p) in NONFH patients versus healthy controls. Likewise, Wei et al.[ screened out 9 up-regulated and 5 down-regulated miRNAs in NONFH serum compared to control serum. Nevertheless, analysis of miRNA expression in TONFH is seldom reported. A recent research screened out a series of key miRNAs involved in TONFH,[ however, miRNA target prediction, function enrichment analysis, and protein-protein network analysis which are helpful for providing a comprehensive understanding of the mechanisms of key miRNAs in TONFH, are still elusive in that study.

In the present study, to identify promising miRNA biomarkers and provide a comprehensive understanding of the mechanisms of promising miRNAs in traumatic osteonecrosis development, we conducted a bioinformatics analysis using the miRNA expression profile dataset from Gene Expression Omnibus. We examined differentially expressed (DE) miRNAs in the blood of TONFH patients and controls with rather strict thresholds. We then performed target prediction, function enrichment analysis, and protein-protein interaction (PPI) network analysis based on DE miRNAs.

Materials and methods

MiRNA expression dataset

The GSE89587[ miRNA expression dataset, containing raw and normalized data, was downloaded from Gene Expression Omnibus database. This dataset consists of 10 samples from TONFH patients and 10 samples from patients who lacked ONFH after trauma. The GPL21439 (miRCURY LNA miRNA Array, 7th generation-hsa) was used as microarray platform. The normalized data was undertaken for further analysis in our study. Approval from a local ethics committee was not required because the GSE89587 database is publicly available.

Identification of DE miRNA

The Linear Models for Microarray Data[ package based on R (version 3.5.1) software and bioconductor[ was undertaken to screen out DE miRNAs between the samples of TONFH and healing patients. Adjusted P-value, also known as false discovery rate (FDR), was analyzed using Banjamini and Hochberg's method.[ DE miRNAs should meet the thresholds of |log2 fold change (FC)|≥ 2.37 and adjusted P-value < .05. The heatmap and volcano plot of DE miRNAs were generated by pheatmap package and ggplot2 package in R software, respectively.

Prediction of target genes

MiRTarBase (https://mirtarbase.mbc.nctu.edu.tw/php/)[ is a freely open and experimentally verified miRNA-target interactions database. MiRDB (https://www.mirdb.org/)[ is a publicly accessible online analysis tool for target prediction of miRNAs as well. Target genes for miRNA were determined in a combination of the 2 databases. First, only the top 3 most up-regulated and down-regulated miRNAs were selected for prediction.[ Second, we only retained validated miRNAs targets with strong evidence within miRTarBase. Third, predicted targets with score > 85 were kept in miRDB.[ We set these rather strict thresholds in an attempt to filter out false positives.

Functional enrichment analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID 6.8, https://david.ncifcrf.gov/) makes a group of functional interpretation tools available to obtain signaling pathways and biological function from lists of genes. To further explore biological process (BP), molecular function (MF) and cellular component (CC) involved in the molecular mechanisms of TONFH, we conducted functional enrichment analysis for these predicted targets, including Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis.

Construction of PPnetwork

The Search Tool for the Retrieval of Interacting Genes (version 10.5, https://www.string-db.org/)[ is a freely online database and search tool for the interactions construction of predicted genes or proteins. Only interactions with a combined score more than 0.4 were mapped. And then, in order to obtain the list of hub proteins, the degree of connectivity of hub genes was calculated and visualized by Cytoscape (version 3.6.1).[ the hub genes were defined as whose node degree greater or equal to 10.

Construction of predicting model

Least absolute shrinkage and selection operator (Lasso) regression analysis is commonly used to determine the correlation between binary variables and the expression level of each gene. In this study, it was applied to select the subset of DE miRNAs that most strongly correlated with the status of femoral head necrosis. In R platform, Lasso regression method in glmnet[ package is employed to calculate the coefficient of each selected miRNA. Then, the predicting model was constructed with respect to the status of femoral head necrosis. The predicting model formula was presented as follows: predicting score = Σiωiχi, in which ωi and χi represent the coefficient and expression level of each selected miRNAs, respectively. The sensitivity and specificity of the predicting model was depicted by a receiver operating characteristic (ROC) curve using the timeROC[ package.

Results

Identification of DE miRNAs

Based on the cut-off criteria, 3 up-regulated and 23 down-regulated miRNAs were identified in TONFH serum versus control serum (Table 1). According to FC, miR-181a-5p (FDR = 4.67E-02), miR-140-5p (FDR = 4.78E-02), and miR-374c-5p (FDR = 4.83E-02) were the top 3 most up-regulated miRNAs in TONFH serum relative to control serum. And miR-4511 (FDR = 1.92E-02), miR-4711-3p (FDR = 2.87E-02), and miR-4479 (FDR = 4.40E-02) were the top 3 most down-regulated ones. The list of these DE miRNAs was then visualized in a heatmap (Fig. 1), in which DE miRNAs were not only ranked by their normalized FCs but also their degrees of similarity. A volcano plot of up-regulated and down-regulated DE miRNAs was shown in Figure 2.

Table 1

Differentially expressed microRNAs in serum from traumatic osteonecrosis of the femoral head patients and controls.

Figure 1

The heatmap of differentially expressed microRNA expression value in traumatic osteonecrosis of the femoral head patients and controls.

Figure 2

Volcano plot of differentially expressed microRNAs. The red dot represents upregulated microRNA and the green dot represents downregulated microRNA.

Targets predication and GO enrichment analysis

We screened out 77 target genes in prediction using the combination of 2 databases, miRTarBase[ and miRDB.[ GO enrichment analysis was conducted on these target genes in 3 categories, including BP, CC and MF. As shown in Figure 3, a total of 63 GO terms were significantly enriched, which mainly involved chromatin organization, chromosome organization, chromatin modification, transcription, and regulation of transcription in the BP category; nucleoplasm part, intracellular organelle lumen, organelle lumen, nuclear lumen, and chromatin in the CC category; and transcription regulator activity, transcription activator activity, ubiquitin-protein ligase activity, chromatin binding, and transcription factor binding in the MF category.

Figure 3

Gene Ontology enrichment analysis of target genes in the categories of biological process, cellular component, and molecular function.

KEGG enrichment analysis

KEGG enrichment analysis was subsequently performed to further explore the enriched pathways of target genes. Interestingly, ubiquitin mediated proteolysis 5 (P = 4.76E-03) was the only enriched pathway, which contained 5 genes, UBE2E3, UBE2K, UBA6, SIAH1, and UBE2B.

Construction of PPI network

To aid in the understanding of the interactions of target genes, PPI network was constructed by the Search Tool For The Retrieval Of Interacting Genes database. As displayed in Figure 4, RBBP7, ACLY, KAT2B, UBE2B, and UBE2K were identified as hub genes. According to the results of GO enrichment analysis, these genes were mainly involved in small conjugating protein ligase activity, ubiquitin-protein ligase activity, and chromatin remodeling complex.

Figure 4

Protein-protein interaction network of target genes in traumatic osteonecrosis of the femoral head constructed from STRING database. The red one represents a hub gene with node degree greater or equal to 10. STRING = search tool for the retrieval of interacting genes.

Using the Lasso regression analysis, an predicting model was established consisting of 5 DE miRNAs (miRNA-3191-5p, miRNA-4511, miRNA-4711-3p, miRNA-4783, and miRNA-122-5p) in terms of the status of femoral head necrosis. The predicting model was constructed as follows: predicting score= (-1.16 × expression level of miRNA-3191-5p) + (-0.76 × expression level of miRNA-4511) + (-0.82 × expression level of miRNA-4711-3p) + (-83.04 × expression level of miRNA-4783) + (-7.02 × expression level of miRNA-122-5p). Utilizing the ROC and the timeROC analysis, we found that the area under the curve (AUC) of the predicting model was 0.98, which means that the power of predicting femoral head necrosis is 98%. (Fig. 5)

Figure 5

The receiver operating characteristic (ROC) curve of predicting model.

Discussion

TONFH is 1 of major complications following surgical treatment of femoral neck fractures. It was reported that the rate of TONFH was 15% to 30% after femur neck fracture.[ Age, timing of surgery, accompanying diseases and presence of displacement are all clinical risk factors affecting the rate of TONFH after femoral neck fracture.[ Biomechanical and biological factors affecting and accelerating fracture healing are also popular topics in research.[ However, the role of miRNA in the development of TONFH after femoral neck fracture remains unclear.

High-throughput genome sequencing, and more recently, RNA-sequencing, as well as bioinformatics technologies have dramatically expanded our understanding of the roles of miRNAs and the molecular mechanisms underlying TONFH. But their functions in TONFH are hindered by the limited number of miRNA expression profiles. We herein, performed a comprehensive bioinformatics analysis based on a miRNA expression profile of TONFH. The results demonstrated that 26 DE miRNAs were identified in TONFH patients compared with patients who lacked ONFH after trauma.

With the help of Lasso regression analysis, we established an predicting model consisting of 5 DE miRNAs (miRNA-3191-5p, miRNA-4511, miRNA-4711-3p, miRNA-4783, and miRNA-122-5p) in terms of the status of femoral head necrosis. Utilizing the ROC and the timeROC analysis, we found that the AUC of the predicting model was 0.98, which means that the power of predicting femoral head necrosis is 98%. This predicting model might help discriminating high-risk patients developing TONFH after femur neck fracture.

Our results not only constructed a predicting model with miRNAs, but also identified DE miRNAs which might serve as new potential molecular treatment targets of TONFH. ONFH is characterized by collapse of the hip cartilage and femoral head, resulting in the secondary osteoarthritis (OA) at the end stage. The role of miR-181a-5p on epithelial-to-mesenchymal transition has been described in OA.[ Okuhara et al[ indicated that miR-181a-5p was highly expressed in the serum collected from OA patients compared with that from controls. Similarly, Nakamura et al[ found that the expression level of miR-181a-5p was improved in trauma-induced knee OA cartilage compared to control cartilage. Histologically, the loss of articular cartilage amongs 1 of the main characteristics of OA. Several studies have underlined the key role of up-regulated miR-181a in the damage of cartilage and bone in OA.[ The destruction of hip cartilage is also essential for the deterioration of ONFH.[ These data suggested that the upregulation of miR-181a-5p might promote osteonecrosis via the damage the hip cartilage in post-traumatic patients.

MiR-140 is of great importance in cartilage development and homeostasis. It has been found to be highly expressed in the cartilage of developing chickens, zebrafish and mice.[ MiR-140 was mostly studied in cancer and OA. In knee OA, the expression of miR-140-5p was decreased in serum and artilage.[ Likewise, the significant reduced expression level of miR-140-5p was found in NONFH from previous literatures.[ In contrast, the serous expression of miR-140-5p was elevated in TONFH in the present study. In glucocorticoids induced ONFH, glucocorticoids could activate a diverse range of signaling pathways in bone that resulted in the up-regulation of E3 ubiquitin ligases.[ Moreover, Nakamura et al[ found that miR-140 was located in 1 intron of ubiquitin ligases (E3 s). Although the role of miR-140 in the effects of glucocorticoids in bone remains poorly understood, this evidence implicates the pivotal roles of miR-140 and ubiquitin proteasome system (UPS) in the development of osteonecrosis. Further researches were still needed to verify the potential role of miR-140-5p in TONFH. Additionally, there is a lack of data regarding whether miR-374c-5p has a strong connection with TONFH; it may be expected to become a promising biomarker in TONFH in the future studies.

UPS was also emphasized in the current study by predicting total 77 target genes including UBE2E3, UBE2K, UBA6, UBE2B, and so on. Afterwards, GO enrichment analysis indicated that these genes were mainly enriched in E3 activity, and ubiquitin (Ub)-dependent protein catabolic process. Moreover, the enriched pathway was related to ubiquitin mediated proteolysis 5, which contained 5 genes, UBE2E3, UBE2K, UBA6, SIAH1, and UBE2B. Meanwhile, UBE2K and UBE2B were regarded as hub genes in the PPI network.

In the last few years, increasing evidence has revealed the mysterious role of UPS in mediating bone metabolism via the induction of osteoblast differentiation and bone formation, as well as the inhibition of osteoclast formation and bone resorption.[ Ub is a highly conserved protein that consists of 76 amino acids. Ubiquitination of a protein requires the concerted action of 3 classes of proteins, a single Ub-activating enzyme (E1), Ub-conjugating enzymes (E2 s), and E3 s.[ Previous observations have revealed that E2 s lie at the heart of ubiquitylation by regulating Ub chain formation on several levels.[ UBE2 is involved in various processes, such as DNA repair, mutagenesis, and cell proliferation.[ Amongst E2 s, UBE2B, also called RAD6, was 1 of the most studied enzymes so far. Recent study has indicated that UBE2B gene was associated with tumorigenesis in breast cancer through interaction with tumor suppressor p53.[ UBE2K, also named Huntingtin-interacting protein 2 (HIP2), E2-25K, UBCH1, is ubiquitously expressed with the highest expression level in brain.[ In Huntington's disease, UBE2K was involved in aggregate formation of expanded polyglutamine proteins and polyglutamine-induced cell death.[ In Alzheimer's disease, UBE2K was revealed as a mediator of Aβ neurotoxicity through inhibition of proteasome activity and intraneuronal accumulation of ubiquitin conjugates.[ To our best of knowledge, studies of UBE2B and UBE2K focused on TONFH have not been published. Thus, further studies were needed to verify the potential roles of Ubs, especially UBE2B and UBE2K in the development of TONFH.

There are still some limitations in this work. First, the predicting model was constructed on a small sample size, therefore the false positive rate of this model might be high. Second, miRNAs results are susceptible to bioinformatics parameters that may vary among clinical parameters, thus it must be validated prospectively for further confirmation. Despite of these drawbacks, our study first established a predicting model with good power to discriminate high-risk patients who might develop TONFH after femur neck fracture. Also, our results identified DE miRNAs which might serve as new potential molecular treatment targets of TONFH.

Conclusions

A predicting model consisting of 5 miRNAs may help discriminating high-risk patients who might develop TONFH after femur neck fracture. Among DE miRNAs, MiR-181a-5p and miR-140-5p may contribute to osteonecrosis development after femur neck fracture via UPS.

Author contributions

Conceptualization: Biaofang Wei.

Methodology: Ning Chen and Bolian Xiao.

Software: Bolian Xiao and Shiying Wang.

Supervision: Biaofang Wei.

Validation: Ning Chen.

Writing – original draft: Ning Chen.