Skip to main content

Bioinformatics analysis of miRNA and mRNA expression profiles to reveal the key miRNAs and genes in osteoarthritis

Abstract

Background

Osteoarthritis (OA) is a chronic degenerative joint disease and the most frequent type of arthritis. This study aimed to identify the key miRNAs and genes associated with OA progression.

Methods

The GSE105027 (microRNA [miRNA/miR] expression profile; 12 OA samples and 12 normal samples) and GSE48556 (messenger RNA [mRNA] expression profile; 106 OA samples and 33 normal samples) datasets were selected from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) and miRNAs (DEMs) were analyzed using the limma and ROCR packages in R, respectively. The target genes that negatively correlated with the DEMs were predicted, followed by functional enrichment analysis and construction of the miRNA-gene and miRNA-transcription factor (TF)-gene regulatory networks. Additionally, key miRNAs and genes were screened, and their expression levels were verified by real-time quantitative reverse transcription polymerase chain reaction (qRT-PCR).

Results

A total of 1696 DEGs (739 upregulated and 957 downregulated) and 108 DEMs (56 upregulated and 52 downregulated) were identified in the OA samples. Furthermore, 56 target genes that negatively correlated with the DEMs were predicted and found to be enriched in three functional terms (e.g., positive regulation of intracellular protein transport) and three pathways (e.g., human cytomegalovirus infection). In addition, three key miRNAs (miR-98-5p, miR-7-5p, and miR-182-5p) and six key genes (murine double minute 2, MDM2; glycogen synthase kinase 3-beta, GSK3B; transmembrane P24-trafficking protein 10, TMED10; DDB1 and CUL4-associated factor 12, DCAF12; caspase 3, CASP3; and ring finger protein 44, RNF44) were screened, among which the miR-7-5p → TMED10/DCAF12, miR-98-5p → CASP3/RNF44, and miR-182-5p → GSK3B pairs were observed in the regulatory network. Moreover, the expression levels of TMED10, miR-7-5p, CASP3, miR-98-5p, GSK3B, and miR-182-5p showed a negative correlation with qRT-PCR verification.

Conclusion

MiR-98-5p, miR-7-5p, miR-182-5p, MDM2, GSK3B, TMED10, DCAF12, CASP3, and RNF44 may play critical roles in OA progression.

Background

As a chronic degenerative joint disease, osteoarthritis (OA) is induced by the destruction of articular cartilage and bone tissue [1]. The most common symptoms of OA include joint pain, joint swelling, stiffness, and limited motion range [2]. OA, which usually occurs in the lower back, fingers, hips, neck, and knees, is mainly caused by abnormal development of the affected limb or joint, joint injury, and genetic factors [3, 4]. At present, the therapeutic options for OA include reduction of joint stress, exercise, pain medications, and support groups [5]. OA is the most frequent type of arthritis and affects approximately 237 million people globally [6]. Therefore, exploring the mechanisms of OA is necessary in order to improve its diagnosis and therapeutic options.

MicroRNAs (miRNAs/miRs) are non-coding RNAs with lengths of 21–24 nucleotides that play crucial roles in the progression of OA by regulating the expression levels of their specific target genes [7]. Studies have demonstrated that miRNAs are useful for the diagnosis and treatment of OA, and the interactions between miRNA and their targets are implicated in the regulation of gene expression and homeostatic pathways in OA [8, 9]. MiR-140 is downregulated in the cartilage of OA patients, and transfection with double-stranded miR-140 may result in decreased matrix metalloproteinase 13 (MMP13) and increased type II procollagen (COL2A1) [10, 11]. MiRNA-155 and miRNA-146a expression levels are upregulated in the later stages of OA, indicating that they may be implicated in the development and progression of OA [12]. By mediating SMAD family member 3 expression, increased miR-16-5p may promote the progression of OA [13]. Dysregulated miR-9, miR-146, and miR-98 function in regulating tumor necrosis factor alpha (TNF-α) production in primary chondrocytes, and miR-9 can suppress the secretion of MMP13 in the chondrocytes of OA patients [14]. Although previous studies have found that these miRNAs and genes are involved in the progression of OA, the pathogenesis of OA is not fully understood.

In this study, the gene and miRNA expression profiles of OA were downloaded. Then, differential expression analysis was conducted to identify differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs). Furthermore, miRNA-gene prediction, functional enrichment analysis, and miRNA-transcription factor (TF)-gene regulatory network analysis were performed, and the key miRNAs and genes in OA were selected and verified. The key miRNAs and genes identified in this study may be candidate diagnostic and therapeutic targets for OA patients.

Methods

Data source

The OA gene expression and miRNA expression profiles were searched using the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) based on the following selection criteria: (1) the dataset included blood or serum samples from OA patients and normal controls and (2) the samples did not undergo any specific treatment, such as drug treatment. Finally, the GSE105027 (miRNA expression profile, platform: GPL21575 Agilent-070156 Human_miRNA_V21.0_Microarray 046064 Feature Number version; 12 serum samples from OA patients and 12 serum samples from normal controls) and GSE48556 (mRNA expression profile, platform: GPL6947 Illumina HumanHT-12 V3.0 expression beadchip; 106 blood samples from OA patients and 33 blood samples from normal controls) datasets were selected.

Differential expression analysis

The limma package in R (version 3.36.3; http://www.bioconductor.org/packages/release/bioc/html/limma.html) [15] was used for normalization processing. The median value of the probes corresponding to the same gene was taken as the expression value of the gene, and the average value of the probes corresponding to the same miRNA was obtained as the expression value of the miRNA. The P value of mRNA differential analysis was corrected for the false discovery rate (FDR), and genes with a FDR < 0.05 were regarded as DEGs. The ROCR package in R (version 1.0-7; https://cran.r-project.org/web/packages/ROCR/index.html) [16] was used to perform receiver operating characteristic (ROC) analysis of miRNAs in OA and normal samples. The area under the ROC curve (AUC) represented the ability of miRNAs to distinguish OA samples from normal samples, and a P value < 0.05 and |AUC-0.5| > 0.3 were set as the criteria for screening DEMs.

Prediction of miRNA-gene pairs

The target genes of the DEMs were searched using the DIANA-TarBase (version 7.0; http://diana.imis.athena-innovation.gr/DianaTools) [17] and miTarBase (version 7.0; http://mirtarbase.mbc.nctu.edu.tw/php/index.php) [18] databases. Target genes appearing in both databases were selected as the target genes of DEMs, and the target genes that negatively correlated with DEMs were further screened [19].

Enrichment analysis

The Gene Ontology (GO) database is valuable for predicting the biological process (BP), molecular function (MF), and cellular component (CC) functions of gene products [20]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database is typically used to predict pathways that involve certain genes [21]. Using the clusterProfiler package (version 3.8.1; http://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) [22], the target genes that negatively correlated with DEMs underwent GO_BP and KEGG enrichment analyses. A FDR < 0.05 was defined as the significant threshold.

Construction of the miRNA-gene network

For the DEGs that negatively correlated with DEMs, the protein-protein interaction (PPI) pairs were extracted from the STRING database (version 10.5; http://string-db.org/) [23], with a STRING score threshold of > 400. Along with the DEM-DEG pairs with negative correlations, the miRNA-gene network was constructed using the Cytoscape software (http://www.cytoscape.org) [24].

Construction of the miRNA-TF-gene regulatory network

Using the chromatin immunoprecipitation sequencing data of the 153 TFs in the Cistrome Data Browser (http://cistrome.org/db) [25], the TFs with a FDR < 0.05 in the differential expression analysis were screened out. Then, bedtools (version 2.17.0; https://bedtools.readthedocs.io/en/latest/) [26] were used to screen the TFs that could target the genes in the miRNA-gene network. Moreover, the miRNA-TF-gene regulatory network was built using the Cytoscape software [24].

Identification of key miRNAs and genes

The top three miRNAs targeting the highest number of genes in the miRNA-TF-gene regulatory network were considered as the key miRNAs. Regarding protein-coding genes in the miRNA-TF-gene regulatory network, searches were conducted in the National Center for Biotechnology Information PubMed database using gene name and “osteoarthritis” as keywords. Genes with search results were considered as key genes. In addition, the expression levels of key genes and the ROC curves of the key miRNAs are analyzed.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR) verification

Based on the results of the above analysis, the expression levels of the three key miRNAs and their target genes (key genes) were verified in six OA blood samples and six normal blood samples (controls). Total RNA was extracted using TRIzol reagent (No. 9109; Takara Bio, Shiga, Japan) according to the manufacturer’s instructions. The concentration and quality were determined using a microplate reader (Infinite M100 PRO; Tecan, Männedorf, Switzerland). RT of total RNA was performed using PrimeScript™RT Master Mix (No. RR036A; Takara Bio), followed by qRT-PCR using Power SYBR Green PCR Master Mix (No. A25742; Thermo Fisher Scientific, Waltham, MA, USA) under the following thermal cycling conditions: 50 °C for 2 min, 95 °C for 10 min, followed by 40 cycles at 95 °C for 10 s and 60 °C for 30 s. The melting curve was analyzed from 60 to 95 °C at an incremental rate of 0.5 °C/10 s. The primer sequences are shown in Supplementary Table 1. Finally, the relative expression levels of miRNAs and target genes were calculated according to the 2−ΔΔCt method [27]. All data are presented as means ± standard deviations. Statistical analysis and graphing were performed using the GraphPad Prism 5 software (GraphPad Software, San Diego, CA, USA). P < 0.05 and P < 0.01 were used to identify results with statistically significant differences and highly significant differences, respectively.

Ethics and consent

This study was approved by the medical ethics committee of Shanghai Fifth People’s Hospital affiliated with Fudan University, and informed consent was obtained from all patients prior to participation. All microarray datasets were published and downloaded from the GEO database; thus, we confirmed that all necessary ethical approvals and informed consent were obtained.

Results

Differential expression analysis

A total of 1696 DEGs (739 upregulated and 957 downregulated) were identified between OA and normal samples (Fig. 1a). Relative to the normal samples, 108 DEMs (56 upregulated and 52 downregulated) were screened in the OA samples (Fig. 1b).

Fig. 1
figure1

Clustering heatmaps. a Heatmap for the differentially expressed genes (DEGs). b Heatmap for the differentially expressed miRNAs (DEMs). In the sample strip, red and green represent osteoarthritis (OA) and control samples, respectively

Prediction of miRNA-gene pairs

The 108 DEMs had 5856 (involving 11,052 miRNA-gene pairs) and 7194 (involving 17,598 miRNA-gene pairs) target genes in the DIANA-TarBase and miTarBase databases, respectively. A total of 755 target genes (involving 846 miRNA-gene pairs) appeared in both databases and thus were selected as the target genes of the DEMs (Fig. 2a). Among the 846 miRNA-gene pairs, 12 miRNAs (4 upregulated and 8 downregulated) were negatively correlated with 56 target genes (including 43 upregulated and 13 downregulated genes) (Fig. 2b).

Fig. 2
figure2

Venn diagram and negatively correlated miRNA-gene pairs. a Overlap of the target genes in the DIANA-TarBase (red) and miTarBase (green) databases. b The 56 target genes that negatively correlated with 12 miRNAs (red circles, blue circles, and inverted triangles represent upregulated genes, downregulated genes, and miRNAs, respectively). A t shape represents the miRNA-gene regulatory relationship

Enrichment analysis

For the 56 target genes that were negatively correlated with 12 DEMs, GO_BP and KEGG enrichment analyses were carried out. GO_BP enrichment analysis showed that three terms (including positive regulation of intracellular protein transport, positive regulation of intracellular transport, and regulation of intracellular protein transport) were enriched (Fig. 3a). Meanwhile, three KEGG pathways were enriched, including human cytomegalovirus infection, thyroid hormone signaling pathway, and ubiquitin-mediated proteolysis (Fig. 3b).

Fig. 3
figure3

Enrichment analysis results. a The Gene Ontology (GO)_biological process (BP) terms enriched for the 56 target genes. b The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched for the 56 target genes

Construction of the miRNA-gene network

After the 56 target genes were mapped into the STRING database, 21 target genes were found to be involved in the PPI pairs. In addition, the 21 target genes were targeted by 7 DEMs. Finally, the miRNA-gene network was visualized, which involved 28 nodes (21 genes and 7 miRNAs) and 39 edges (21 miRNA-gene regulatory pairs and 18 PPI pairs) (Fig. 4a).

Fig. 4
figure4

Regulatory networks. a The miRNA-gene network. b The miRNA-transcription factor (TF)-gene regulatory network. Circles, inverted triangles, and diamonds represent genes, miRNAs, and TFs, respectively. Red and blue represent upregulated and downregulated genes, respectively. T shapes, lines, and arrows represent miRNA-gene regulatory relationships, protein-protein interactions, and TF-gene regulatory relationships, respectively

Construction of the miRNA-TF-gene regulatory network

A total of 17 TFs were found to be differentially expressed between the OA and normal samples, among which 16 TFs could target the genes in the miRNA-gene network. In the miRNA-TF-gene regulatory network, there were 44 nodes (21 genes, 7 miRNAs, and 16 TFs) and 280 edges (21 miRNA-gene regulatory pairs, 18 PPI pairs, and 241 TF-gene regulatory pairs) (Fig. 4b).

Identification of key miRNAs and genes

MiR-98-5p, miR-7-5p, and miR-182-5p were the top three miRNAs that targeted the highest number of genes in the miRNA-TF-gene regulatory network and thus were considered as the key miRNAs. Six key genes (including murine double minute 2, MDM2; glycogen synthase kinase 3-beta, GSK3B; transmembrane P24-trafficking protein 10, TMED10; DDB1 and CUL4-associated factor 12, DCAF12; caspase 3, CASP3; and ring finger protein 44, RNF44) had search results in the PubMed database. The ROC curves of the three key miRNAs and the expression levels of the six key genes are shown in Figs. 5 and 6, respectively. In particular, the miR-7-5p → TMED10/DCAF12, miR-98-5p → CASP3/RNF44, and miR-182-5p → GSK3B regulatory pairs as well as the MDM2─GSK3B/CASP3 PPI pairs were present in the miRNA-TF-gene regulatory network.

Fig. 5
figure5

Receiver operating characteristic (ROC) curves of the three key miRNAs. a The ROC curve of miR-98-5p. b The ROC curve of miR-7-5p. c The ROC curve of miR-182-5p. AUC, area under the ROC curve

Fig. 6
figure6

Expression levels of the six key genes (murine double minute 2, MDM2; glycogen synthase kinase 3-beta, GSK3B; transmembrane P24-trafficking protein 10, TMED10; DDB1 and CUL4-associated factor 12, DCAF12; caspase 3, CASP3; and ring finger protein 44, RNF44)

Expression levels of the key miRNAs and genes

The expression levels of the three key miRNAs (miR-98-5p, miR-7-5p, and miR-182-5p) and three key target genes (TMED10, CASP3, and GSK3B) were verified using qRT-PCR. As shown in Supplementary Figure 1, the OA samples displayed a significant decrease in the expression levels of miR-7-5p and miR-98-5p and a significant increase in the expression level of miR-182-5p compared to those of the control samples. In addition, the expression levels of TMED10, CASP3, and GSK3B in the OA samples changed as expected. More specifically, the TMED10 (the target gene of miR-7-5p) and CASP3 (the target gene of miR-98-5p) expression levels were significantly upregulated while that of GSK3B (the target gene of miR-182-5p) was significantly downregulated in the OA samples compared to the controls.

Discussion

In this study, 1696 DEGs (739 upregulated and 957 downregulated) and 108 DEMs (56 upregulated and 52 downregulated) were screened in OA samples relative to normal samples. Among the predicted miRNA-gene pairs, 56 target genes were negatively correlated with 12 DEMs. For the 56 target genes, three GO_BP terms (such as positive regulation of intracellular protein transport) and three KEGG pathways (such as human cytomegalovirus infection) were enriched. In addition, the miRNA-gene and miRNA-TF-gene regulatory networks were constructed. Moreover, three key miRNAs (miR-98-5p, miR-7-5p, and miR-182-5p) and six key genes (MDM2, GSK3B, TMED10, DCAF12, CASP3, and RNF44) were screened, among which the miR-7-5p → TMED10/DCAF12, miR-98-5p → CASP3/RNF44, and miR-182-5p → GSK3B regulatory pairs were found to be present in the regulatory network. In particular, the expression levels of TMED10 and miR-7-5p, CASP3, miR-98-5p, GSK3B, and miR-182-5p showed a negative correlation with qRT-PCR verification.

MiR-7-5p and miR-214-5p expression levels are elevated in rheumatoid arthritis (RA) with interstitial lung disease (ILD), which may serve as diagnostic markers for RA-associated ILD [28]. TMED10 expression is critical for transforming growth factor β signaling [29], which influences cartilage integrity as well as the prevention and repair of cartilage damage [30]. DCAF12 belongs to the WD40-motif repeat protein family, which has evolutionary conservation and exerts multiple roles in signal transduction [31]. The miR-7-5p → TMED10/DCAF12 regulatory pairs existed in the miRNA-TF-gene regulatory network, indicating that miR-7-5p may act by regulating TMED10 and DCAF12 in OA.

MiR-98 expression is decreased in OA chondrocytes, and its inhibition can result in chondrocyte apoptosis; therefore, miR-98 may be applied in targeted OA therapy [32, 33]. By enhancing extracellular signal-regulated kinase expression and inhibiting CASP3, FAS, and Fas ligand expression, recombinant human lactoferrin can reverse dexamethasone-induced chondrocyte apoptosis in OA [34]. Upregulated Piezo1 affects the apoptotic rate of OA chondrocytes via a CASP12-dependent pathway [35]. Through the histone deacetylase 9 protein inhibitor of the activated STAT Y-RNF4 axis, the post-translational modification of Nkx3.2 acts to mediate chondrocyte hypertrophy during the course of skeletal development [36]. RNF4 is targeted by Epstein-Barr virus miR-BHRF1-1, and its reconstitution is related to decreased viral proteins and damaged virus release [37]. The MiR-98-5p → CASP3/RNF44 regulatory pairs were involved in the miRNA-TF-gene regulatory network, suggesting that miR-98-5p may affect the pathogenesis of OA by targeting CASP3 and RNF44.

MiR-182 positively mediates T helper cell proliferation by suppressing forkhead box O1 expression, and its inhibition prevents T helper cells from inducing arthritis [38]. GSK3B maintains the chondrocytic phenotype and cartilage extracellular matrix by suppressing β-catenin, and its inhibition promotes the differentiation of OA chondrocytes [39]. The MiR-182-5p → GSK3B regulatory pair was found in the miRNA-TF-gene regulatory network; therefore, miR-182-5p may also be implicated in the development of OA by mediating GSK3B. MDM2 expression in the fibroblast-like synoviocytes of RA patients was upregulated compared to that of OA patients, which may help to reduce apoptosis of lining tissues in RA [40]. MDM2 interacted with GSK3B and CASP3 in the miRNA-TF-gene regulatory network; thus, MDM2 may function in the OA mechanisms via interaction with GSK3B and CASP3.

Despite these findings, this study has some limitations. The lack of human disease tissues that reflect disease characteristics limits the microarray study of human diseases, including OA. Considering the roles of peripheral blood in immune response, metabolism, and intercellular communication, especially the simplicity of sample collection, it has been considered as a notable tissue for biomarker detection. Therefore, a microarray dataset of blood samples was selected and analyzed. However, blood samples could be affected by many factors, and the data from tissue samples might provide more reliable results. This is one limitation of the present study. In addition, we only confirmed the mRNA expression levels of the key genes and miRNAs. MiRNA-mRNA interactions should be further investigated using a dual-luciferase reporter assay. Moreover, the protein expression levels of key target genes should be further confirmed by western blot analysis. Finally, the roles of these key genes and miRNAs in the development and progression of OA should be investigated using functional experiments.

In conclusion, 1696 DEGs and 108 DEMs were identified between OA and normal samples. Additionally, miR-98-5p, miR-7-5p, miR-182-5p, MDM2, GSK3B, TMED10, DCAF12, CASP3, and RNF44 may play important roles in the pathogenesis of OA. However, the functions of these key miRNAs and genes in OA need to be validated by future experiments.

Availability of data and materials

The datasets used and analyzed during the current study are available from the corresponding author upon reasonable request.

Abbreviations

BP:

Biological process

CC:

Cellular component

DEGs:

Differentially expressed genes

DEMs:

Differentially expressed miRNAs

GEO:

Gene Expression Omnibus

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MF:

Molecular function

OA:

Osteoarthritis

RA:

Rheumatoid arthritis

TF:

Transcription factor

References

  1. 1.

    Jie S, Di C. Recent progress in osteoarthritis research. J Am Acad Orthop Surg. 2014;22(7):467–8.

    Article  Google Scholar 

  2. 2.

    Owens C, Conaghan PG. Improving joint pain and function in osteoarthritis. Practitioner. 2016;260(1799):17.

    PubMed  Google Scholar 

  3. 3.

    Li G, Yin J, Gao J, Cheng TS, Pavlos NJ, Zhang C, et al. Subchondral bone in osteoarthritis: insight into risk factors and microstructural changes. Arthritis Res Ther. 2013;15(6):223.

    CAS  Article  Google Scholar 

  4. 4.

    Heidari B. Knee osteoarthritis prevalence, risk factors, pathogenesis and features: part I. Caspian J Intern Med. 2011;2(2):205.

    PubMed  PubMed Central  Google Scholar 

  5. 5.

    Mcalindon TE, Bannuru RR, Sullivan MC, Arden NK, Berenbaum F, Biermazeinstra SM, et al. OARSI guidelines for the non-surgical management of knee osteoarthritis. Osteoarthr Cartil. 2014;22(3):363–88.

    CAS  Article  Google Scholar 

  6. 6.

    March L, Smith EUR, Hoy DG, Cross MJ, Sanchez-Riera L, Blyth F, et al. Burden of disability due to musculoskeletal (MSK) disorders. Best Pract Res Clin Rheumatol. 2014;28(3):353–66.

    Article  Google Scholar 

  7. 7.

    Trachana V, Ntoumou E, Anastasopoulou L, Tsezou A. Studying microRNAs in osteoarthritis: critical overview of different analytical approaches. Mech Ageing Dev. 2018;171:15–23. https://pubmed.ncbi.nlm.nih.gov/29496549/.

  8. 8.

    Oliviero A, Della Porta G, Peretti GM, Maffulli N. MicroRNA in osteoarthritis: physiopathology, diagnosis and therapeutic challenge. Br Med Bull. 2019;130(1):137–47.

    CAS  Article  Google Scholar 

  9. 9.

    Giordano L, Porta GD, Peretti GM, Maffulli N. Therapeutic potential of microRNA in tendon injuries. Br Med Bull. 2020;133(1):79–94.

    CAS  Article  Google Scholar 

  10. 10.

    Zhou XH, Wang M, Yan-Hui JI, Jiang Y, Liu Q, Orthopaedics DO, et al. Expression of miRNA-140 in chondrocytes of patients with early osteoarthritis and its function. Acad J Second Military Med Univ. 2014;35(6):610–5.

    CAS  Article  Google Scholar 

  11. 11.

    Si HB, Zeng Y, Zhou ZK, Pei FX, Lu YR, Cheng JQ, et al. Expression of miRNA-140 in chondrocytes and synovial fluid of knee joints in patients with osteoarthritis. Chin Med Sci J. 2016;31(4):207–12.

    Article  Google Scholar 

  12. 12.

    Soyocak A, Kurt H, Ozgen M, Cosan DT, Colak E, Gunes HV. miRNA-146a, miRNA-155 and JNK expression levels in peripheral blood mononuclear cells according to grade of knee osteoarthritis. Gene. 2017;627:207–11.

    CAS  Article  Google Scholar 

  13. 13.

    Lisong L, Jie J, Xianzhe L, Shuhua Y, Shunan Y, Wen Y, et al. MicroRNA-16-5p controls development of osteoarthritis by targeting SMAD3 in chondrocytes. Curr Pharm Des. 2015;21(35):5160–7.

  14. 14.

    Jones S, Watkins GG, Le Good N, Roberts S, Murphy C, Brockbank S, Needham M, et al. The identification of differentially expressed microRNA in osteoarthritic tissue that modulate the production of TNF-alpha and MMP13. Osteoarthr Cartil. 2009;17(4):464–72.

    CAS  Article  Google Scholar 

  15. 15.

    Smyth GK. limma: linear models for microarray data. Bioinformatics & Computational Biology Solutions Using R & Bioconductor; 2011. p. 397–420.

    Google Scholar 

  16. 16.

    Tobias S, Oliver S, Niko B, Thomas L. ROCR: visualizing classifier performance in R. Bioinformatics. 2005;21(20):3940–1.

  17. 17.

    Vlachos IS, Paraskevopoulou MD, Dimitra K, Georgios G, Thanasis V, Ilias K, et al. DIANA-TarBase v7.0: indexing more than half a million experimentally supported miRNA:mRNA interactions. Nucleic Acids Res. 2015;43(Database issue):153–9.

    Article  Google Scholar 

  18. 18.

    Chou CH, Shrestha S, Yang CD, Chang NW, Lin YL, Liao KW, et al. miRTarBase update 2018: a resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 2018;46(D1):D296–302.

  19. 19.

    Fu J, Tang W, Du P, Wang G, Chen W, Li J, et al. Identifying microRNA-mRNA regulatory network in colorectal cancer by a combination of expression profile and bioinformatics analysis. BMC Syst Biol. 2012;6(1):68.

    CAS  Article  Google Scholar 

  20. 20.

    Consortium TGO. Gene Ontology Consortium: going forward. Nucleic Acids Res. 2015;43(Database issue):1049–56.

    Article  Google Scholar 

  21. 21.

    Kanehisa M, Sato Y, Kawashima M, Furumichi M, Tanabe M. KEGG as a reference resource for gene and protein annotation. Nucleic Acids Res. 2016;44(Database issue):D457–D62.

    CAS  Article  Google Scholar 

  22. 22.

    Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    CAS  Article  Google Scholar 

  23. 23.

    Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, et al. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(Database issue):D362–D8.

    CAS  Article  Google Scholar 

  24. 24.

    Kohl M, Wiese S, Warscheid B. Cytoscape: software for visualization and analysis of biological networks. Methods Mol Biol. 2011;696(696):291–303.

    CAS  Article  Google Scholar 

  25. 25.

    Mei S, Qin Q, Wu Q, Sun H, Zheng R, Zang C, et al. Cistrome Data Browser: a data portal for ChIP-Seq and chromatin accessibility data in human and mouse. Nucleic Acids Res. 2017;45(Database issue):D658–D62.

    CAS  Article  Google Scholar 

  26. 26.

    Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010.

  27. 27.

    Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2− ΔΔCT method. Methods. 2001;25(4):402–8.

    CAS  Article  Google Scholar 

  28. 28.

    Oka S, Furukawa H, Shimada K, Hashimoto A, Komiya A, Fukui N, et al. Plasma miRNA expression profiles in rheumatoid arthritis associated interstitial lung disease. BMC Musculoskelet Disord. 2017;18(1).

  29. 29.

    Nakano N, Tsuchiya Y, Kako K, Umezaki K, Sano K, Ikeno S, et al. TMED10 interferes with TGF-β signaling by disrupting TGF-β receptor complex formation. J Biol Chem. 2017;292(10):4099–112.

  30. 30.

    Davidson ENB, Kraan PMVD, Berg WBVD. TGF-β and osteoarthritis. Osteoarthritis\s&\scartilage. 2007;15(6):597–604.

    Google Scholar 

  31. 31.

    Xu C, Min J. Structure and function of WD40 domain proteins. Protein Cell. 2011;2(3):202–14.

    CAS  Article  Google Scholar 

  32. 32.

    Wang J, Chen L, Jin S, Lin J, Zheng H, Zhang H, et al. Altered expression of microRNA-98 in IL-1Î2-induced cartilage degradation and its role in chondrocyte apoptosis. Mol Med Rep. 2017;16(3):3208–16.

    CAS  Article  Google Scholar 

  33. 33.

    Wang GL, Wu YB, Liu JT, Li CY. Upregulation of miR-98 inhibits apoptosis in cartilage cells in osteoarthritis. Genet Test Mol Biomark. 2016;20(11):645.

    CAS  Article  Google Scholar 

  34. 34.

    Tu Y, Xue H, Francis W, Davies AP, Pallister I, Kanamarlapudi V, et al. Lactoferrin inhibits dexamethasone-induced chondrocyte impairment from osteoarthritic cartilage through up-regulation of extracellular signal-regulated kinase 1/2 and suppression of FASL, FAS, and Caspase 3. Biochem Biophys Res Commun. 2013;441(1):249–55.

    CAS  Article  Google Scholar 

  35. 35.

    Li XF, Zhang Z, Chen ZK, Cui ZW, Zhang HN. Piezo1 protein induces the apoptosis of human osteoarthritis-derived chondrocytes by activating caspase-12, the signaling marker of ER stress. Int J Mol Med. 2017;40(3):845–53.

    CAS  Article  Google Scholar 

  36. 36.

    Choi HJ, Kwon S, Kim DW. A post-translational modification cascade employing HDAC9-PIASy-RNF4 axis regulates chondrocyte hypertrophy by modulating Nkx3.2 protein stability. Cell Signal. 2016;28(9):1336–48.

    CAS  Article  Google Scholar 

  37. 37.

    Li J, Callegari S, Masucci MG. The Epstein-Barr virus miR-BHRF1-1 targets RNF4 during productive infection to promote the accumulation of SUMO conjugates and the release of infectious virus. PLoS Pathog. 2017;13(4):e1006338.

    Article  Google Scholar 

  38. 38.

    Stittrich AB, Haftmann C, Sgouroudis E, Fang Z, Rajewsky N, Chang HD, et al. Inhibition of miR-182 decreases OVA-induced arthritis in mice. Ann Rheum Dis. 2011;70(Suppl 2):A43.

  39. 39.

    Minguzzi M, Guidotti S, Platano D, Olivotto E, Facchini A, Flamigni F, et al. GSK3β inhibition induces terminal differentiation and extracellular matrix remodelling in human osteoarthritic articular chondrocytes. Osteoarthritis Cartilage. 2014;22(Suppl):S135–S6.

    Article  Google Scholar 

  40. 40.

    Elliott T, Rong XJ, Derek L, Paul H, Malcolm S, Morand EF, et al. Detection of the p53 regulator murine double-minute protein 2 in rheumatoid arthritis. J Rheumatol. 2005;32(3):424.

    Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was supported by the Scientific Research Project funded by Shanghai Fifth People's Hospital, Fudan University (No. 2020WYZT06).

Author information

Affiliations

Authors

Contributions

Conception and design of the research: PH and LL; acquisition of the data: JW and LL; analysis and interpretation of the data: JG and MW; statistical analysis: TZ and SW; drafting of the manuscript: PH; revision of the manuscript for important intellectual content: MW. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Minghai Wang.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the medical ethics committee of Shanghai Fifth People’s Hospital affiliated with Fudan University. Informed consent was obtained from all patients prior to participation.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Supplementary Table 1.

Primer sequences of key miRNAs and genes used in qRT-PCR.

Additional file 2: Supplementary Figure 1.

Expression levels of three key miRNAs and three key target genes determined by qRT-PCR. MiR-98-5p and miR-7-5p were downregulated while miR-182-5p was upregulated in OA samples compared to the controls (P < 0.01). MDM2 and CASP3 were upregulated while GSK3B was downregulated in OA samples compared to the controls (P < 0.01).

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Huang, Py., Wu, Jg., Gu, J. et al. Bioinformatics analysis of miRNA and mRNA expression profiles to reveal the key miRNAs and genes in osteoarthritis. J Orthop Surg Res 16, 63 (2021). https://doi.org/10.1186/s13018-021-02201-2

Download citation

Keywords

  • Osteoarthritis
  • Differentially expressed genes
  • Differentially expressed miRNAs
  • Transcription factors
  • Regulatory network