Skip to main content

Identification of transcription factors and construction of a novel miRNA regulatory network in primary osteoarthritis by integrated analysis

Abstract

Backgrounds

As osteoarthritis (OA) disease-modifying therapies are not available, novel therapeutic targets need to be discovered and prioritized. Here, we aim to identify miRNA signatures in patients to fully elucidate regulatory mechanism of OA pathogenesis and advance in basic understanding of the genetic etiology of OA.

Methods

Six participants (3 OA and 3 controls) were recruited and serum samples were assayed through RNA sequencing (RNA-seq). And, RNA-seq dataset was analysed to identify genes, pathways and regulatory networks dysregulated in OA. The overlapped differentially expressed microRNAs (DEMs) were further screened in combination with the microarray dataset GSE143514. The expression levels of candidate miRNAs were further validated by quantitative real-time PCR (qRT-PCR) based on the GEO dataset (GSE114007).

Results

Serum samples were sequenced interrogating 382 miRNAs. After screening of independent samples and GEO database, the two comparison datasets shared 19 overlapped candidate micRNAs. Of these, 9 up-regulated DEMs and 10 down-regulated DEMs were detected, respectively. There were 236 target genes for up-regulated DEMs and 400 target genes for those down-regulated DEMs. For up-regulated DEMs, the top 10 hub genes were KRAS, NRAS, CDC42, GDNF, SOS1, PIK3R3, GSK3B, IRS2, GNG12, and PRKCA; for down-regulated DEMs, the top 10 hub genes were NR3C1, PPARGC1A, SUMO1, MEF2C, FOXO3, PPP1CB, MAP2K1, RARA, RHOC, CDC23, and CREB3L2. Mir-584-5p-KRAS, mir-183-5p-NRAS, mir-4435-PIK3R3, and mir-4435-SOS1 were identified as four potential regulatory pathways by integrated analysis.

Conclusions

We have integrated differential expression data to reveal putative genes and detected four potential miRNA-target gene pathways through bioinformatics analysis that represent new mediators of abnormal gene expression and promising therapeutic targets in OA.

Peer Review reports

Background

Osteoarthritis (OA), the most common form of arthritis, is a degenerative joint disease characterized by progressive articular cartilage and a leading cause of disability in the elderly [1, 2]. OA affects more than 300 million people according to the recent Global Burden of Disease study 2017, which is the third most rapidly rising condition associated with disability after diabetes and dementia [1, 3]. It is estimated that nearly 400 million people will suffer from OA by 2030 in China [4]. The etiology of OA is affected by factors including age, female, obesity, occupational exposure to high levels of joint loading activity, smoking status and family history, and the molecular pathophysiology is still poorly understood [4]. Epidemiological studies have reported the heritability of >40% at individual skeletal sites, while genome-wide association studies (GWAS) have revealed that OA is a heterogeneous chronic disease with multiple genetics [5]. As such, the disease is genetically complex and multi-factorial.

MicroRNA (miRNA) are a group of universally present and multifunctional small non-coding RNAs with 22–28 bases encoded by endogenous genes [6, 7], and participate in many biological processes such as cell proliferation, differentiation, apoptosis [8]. Recent breakthroughs have revealed numerous examples of miRNAs involvement in normal development and disease [9]. As one of the epigenetic mechanisms, miRNAs play a significant role in OA pathogenesis, such as cartilage homeostasis, extracellular matrix regulation, bone metabolism, apoptosis and inflammation [10]. And, much attention has recently focused on miRNA because increasing evidence indicates that miRNA affect gene transcription through a number of regulatory processes [10].

Mature miRNAs are structurally stable and essential regulators of many physiological processes. Several studies have identified numerous differentially expressed microRNAs (DEMs) between OA joint tissues and controls [11, 12]. Presumable studies demonstrated the role of miRNAs as potential biomarkers for the assessment of diagnosis and prognosis in OA [13, 14]. Nakamura et al. showed that miR-181a-5p induces articular cartilage degeneration by promoting inflammation, cartilage catabolism and apoptosis/cell death [15]. Ntoumou et al. identified that a three-miRNA signature, hsa-miR-140-3p, hsa-miR-681-3p and hsa-miR-33b-3p in OA patients, which were involved in the regulation of OA metabolism [16]. Presently, medical care is mainly based on alleviating pain symptoms and there is still a lack of effective means for the treatment of disease. It is necessary to fully understand the molecular and cellular pathways of OA and novel breakthrough in the sector of early diagnosis and treatment of OA.

In the present study, we performed miRNA profiling using high-density miRNA-arrays in peripheral blood leukocytes (PBL) and compared the expression profiles of miRNAs between OA patients and controls. Then, overlapped DEMs were screened by downloading the miRNAs expression data from the Gene Expression Omnibus (GEO) database (GSE143514). We created a molecular profile of the OA transcriptome and performed functionl enrichment analyses to elucidate perturbed molecular functions and pathways in OA, including transcription factors (TFs)-DEMs, DEMs-target genes, TFs of target genes, and DEMs-hub genes network. Hub genes were further performed to validate the expression levels based on the GSE114007 dataset. Finally, we developed a high-resolution molecular profile of OA and elucidated a novel miRNA regulatory network of OA pathogenesis.

Materials and methods

Sample collection

A two-stage study was conducted in this study. In the first RNA-sequencing screening stage, blood samples were collected form 3 OA patients and 3 healthy controls (age over 40 years old) in March 2018.

The inclusion criteria for OA were presented as follows: (i) Complete basic personal information; (ii) Conforming to the diagnosis criteria of disease [17]. The normal controls were recruited from the survey that had no signs or symptoms of arthritis or joint diseases. Secondary OA patients such as inflammatory arthritis, rheumatoid, bone fracture and developmental dysplasia were excluded. MicroRNA expression dataset was downloaded from the GEO database as another screening stage under the accession number GSE143514.

Subsequently, in the second stage validation, the expression levels of hub genes were further verified using quantitative reverse transcription polymerase chain reaction (qRT-PCR) based on the GSE114007 dataset.

Library preparation and sequencing

White blood cells are extracted by centrifugation at 1500 g for 20 min with 5 ml whole blood. Following the manufacture’s protocol, we used TRIzol (Invitrogen, Carlsbad, CA, USA) to isolate total RNA.

RNA-seq was operated by Gminix Biotechnology Co., Ltd. (Shanghai, China) using 150 ng of total RNA as input, and the results was evaluated by an Illumina HiSeq 2500 sequencing platform with 10 M reads (Illumina, San Diego, CA, USA). Hisat2 [18] were used to compare DEMs, and feature Counts [19] were adopted to annotate and quantify miRNAs. MiRNAs expression was normalized by principal component analysis (PCA) using R software. Then, using the DESeq2 package in R (http://bioconductor.org/packages/release/bioc/html/DESeq2.html), counts were analyzed and DEMs were screened. Herein, miRNAs with an absolute value fold changes (FC) ≥ 1.2 and P < 0.05 were considered as the significance.

Data Collection for GSE143514 and Identification of candidate DEMs

The GEO database (https://www.ncbi.nlm.nih.gov/gds) is an international public repository that provides high-throughput gene expression data [20]. The GEO dataset GSE143514, including five OA patients and three controls were screened [21]. All samples were detected by the platform HiSeq x-ten platform (Illumina). For data processing, the raw data were downloaded (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE143514). The limma software package in R (https://bioconductor.org/packages/release/bioc/html/limma.html) was employed to screen DEMs for OA patients and healthy controls. The absolute value fold changes (FC) ≥1.2 and P < 0.05 were set as the threshold for identifying DEMs.

Respectively, the volcano maps of DEMs in the internal dataset and GSE143514 database were drawn using GraphPad Prism. Venn diagram was utilized to analyze the overlapped DEMs between the internal dataset and GEO dataset (GSE143514), and these miRNAs were considered as candidate DEMs.

Prediction of Upstream Transcription Factors and Downstream Target Genes of DEMs

FunRich, an independent software tool for the analysis of gene and protein functional enrichment and interaction networks [22], was employed to predict the potential upstream TFs of candidate DEMs. Finally, we presented the top 10 transcription factors. The P value<0.05 was considered statistically significant.

Based on the above primary analysis of candidate DEMs, target genes were synchronously predicted using miRTarBase 7.2 (http://www.targetscan.org/vert_72/), miRDB (http://mirdb.org/), miRwalk (http://mirwalk.umm.uni-heidelberg.de/) and MicroT-CDS (http://diana.imis.athena-innovation.gr/DianaTools/index.php?r=microT_CDS/index). For a relatively robust selection of target genes, the overlapped target genes from the above four databases were detected by Venn diagram analysis.

Functional Analysis of Target Genes

Bioinformatics analysis was used to investigate the involved pathways and target genes for the above miRNAs. Gene Ontology (GO) is a comprehensive gene annotation tool that provides information on gene function from molecular level to organism point, including biological process (BP), molecular function (MF), and cellular component (CC) [23]. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a comprehensive database that broadly interprets biological pathways and identifies the functions of candidate DEMs in biological pathways [24]. In our study, the database for annotation, visualization, and integrated discovery (DAVID) v6.8 (https://david.ncifcrf.gov) online tool was performed for candidate DEMs enrichment analysis. P < 0.05 for GO analysis and P < 0.1 and count>2 for the KEGG pathway were considered for further analysis.

Construction of PPI Network and Screening of Hub Genes

Investigation on the key miRNA and target genes protein-protein interaction (PPI) network is beneficial for the exploration of molecular mechanism of the disease. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) online database (http://stringdb.org/) was employed to analyze the interaction between miRNA and target genes. PPI node pairs with a combined score ≥ 0.4 were selected for further analysis. We further matched up-regulated target genes with down-regulated target genes and down-regulated miRNA with up-regulated genes. Cytoscape 3.7.1 was used to visualize the network and annotate genes within enriched pathways in the network. The hub genes were screened by CytoHubba, a plugin of Cyoscape, according to the degree. The top 30 target genes were selected using the Maximal Clique Centrality (MCC).

Expression Analysis of Hub Genes based on GSE114007

For further verification analysis, the GSE114007 database was introduced to analyze the expression level of hub genes. The dataset was performed on the platform GPL11154 Illumina HiSeq 2000 and GPL18573 Illumina NextSeq 500, including 20 OA knee cartilages and 18 normal cartilage tissues. DEMs between OA and control group were analyzed by Student’s t-test. The hub genes must conform to the following criteria: (i) the relationship of gene regulation is definitely consistent with previous studies; (ii) the P value should be less than 0.05.

Statistical analysis

Body mass index (BMI) was calculated by dividing the weight (kg) by the squared height (m2), and the participants were divided into four categories according to the Working Group on Obesity in China recommended criteria (Underweight: BMI < 18.5 kg/m2; Normal: 18.5 ≤ BMI < 24; Overweight: 24 ≤ BMI < 28; and General obesity: BMI ≥ 28), [25]. A value of P < 0.05 was considered statistically significant. All analysis was performed with SPSS 22.0 (IBM Corp) and GraphPad Prism version 9.0 (GraphPad Software, lnc).

Results

A schematic diagram of the study design is displayed in Fig. 1, and the characteristics of the internal database are shown in Table 1.

Fig. 1
figure 1

Flow chart of constructing the miRNA regulatory network in OA

Table 1 Characteristics of the subjects enrolled for miRNA expression analysis in the study

Identification of DEMs

Six serum samples were sequenced based on the selection criteria of adjusted P < 0.05 and |logFC| > 1.2. Briefly, a total of 382 DEMs were identified from the internal database, including 194 up-regulated miRNAs and 188 down-regulated miRNAs. Similarly, 128 DEMs (47 up-regulated and 81 down-regulated) were identified in the GSE143514 dataset. The volcano maps of the two databases were presented in Fig. 2a and b, respectively.

Fig. 2
figure 2

Identification of differentially expressed miRNAs (DEMs). a DEMs of the GSE143514 database. b DEMs of the internal database. c intersection of the two databases with down-regulated DEMs. d intersection of the two databases with up-regulated DEMs

` After screening of independent samples and GEO database, the two comparison datasets shared 19 overlapped candidate miRNAs. The co-existing DEMs were displayed in Fig. 2c and d via Venn diagrams. Among these DEMs, we detected 9 up-regulated DEMs (hsa-mir-34a-5p, hsa-mir-4777-3p, hsa-mir-129-2-3p, hsa-mir-152-5p, hsa-mir-200a-3p, hsa-mir-6503-5p, hsa-mir-138-5p, hsa-mir-501-5p, and hsa-mir-342-5p) and 10 down-regulated DEMs (hsa-mir-185-5p, hsa-mir-4732-3p, hsa-mir-4732-5p, hsa-mir-3143, hsa-mir-584-5p, hsa-mir-4435, hsa-mir-1246, hsa-mir-151a-3p, hsa-mir-183-5p, and hsa-mir-501-3p).

Prediction of Transcription Factors of DEMs

In the study, we used FunRich software to predict the upstream transcription factors of up-regulated and down-regulated DEMs and presented the top 10 transcription factors in Fig. 3a and b, respectively. For up-regulated DEMs, the top 10 transcription factors include EGR1, SP1, POU2F1, SP4, NFIC, RREB1, TCF3, TFAP4, E2F1, and FOXD3. For down-regulated DEMs, the top 10 transcription factors were EGR1, NKX6–1, SP1, POU2F1, HOXA3, FOXK1, HOXB13, POU6F1, SP4, and NFIC.

Fig. 3
figure 3

Predicted transcription factors of DEMs. a top ten upstream transcription factors of up-regulated DEMs. b top ten upstream transcription factors of down-regulated DEMs

Prediction of Target Genes of DEMs

To identify the miRNA-mRNA regulatory axes involved in OA, we used miRTarBase 7.2, miRDB, miRwalk and MicroT-CDS databases to predict the target genes of DEMs. As shown in Table 2, we explored the target genes of miRNAs based on different databases. The overlapped target genes from the above four databases were screened by Venn diagram analysis and presented in Supplementary Table 1. There were 236 target genes for up-regulated DEMs and 400 ones for down-regulated DEMs. The DEMs-target gene network were plotted to more intuitively explore the relationship between DEMs and target genes (Fig. 4a and b).

Table 2 miRNA-target genes count
Fig. 4
figure 4

Potential target genes of DEMs predicted by databases. a up-regulated miRNA-target gene network. b down-regulated miRNA-target gene network

Enrichment Analysis of Target Genes

We further explored biological functions of those candidate target genes through the DAVID database. The top 10 GO enrichment results of up-regulated and down-regulated target genes were shown in Fig. 5a and b, respectively. In the BP group, the up-regulated genes were significantly enriched in gene expression, regulation of nitrogen compound metabolic process, RNA metabolic process, regulation of gene expression, and regulation of nucleobase-containing compound metabolic process, while the down-regulated genes were enriched in regulation of nitrogen compound metabolic process, positive regulation of metabolic process, regulation of cell communication, and nervous system development, etc. In the CC group, the up-regulated genes were significantly enriched in nucleoplasm, neuron part, ribonucleoprotein complex, intracellular ribonucleoprotein complex, and nucleoplasm part, while the down-regulated genes were mainly involved the neuron part, cell junction, plasma membrane region, neuron projection, and synapse, etc. In the MF group, the up-regulated genes were significantly enriched in organic cyclic compound binding, heterocyclic compound binding, nucleic acid binding, RNA binding, and nucleic acid binding transcription factor activity, while the down-regulated genes were mainly involved in organic cyclic compound binding, heterocyclic compound binding, small molecule binding, enzyme binding, and purine ribonucleotide binding, etc.

Fig. 5
figure 5

GO annotation analysis for the target genes of DEMs in the biological process, cellular component, and molecular function. a GO analysis for up-regulated DEG. b GO analysis for down-regulated DEGs

Subsequently, we conducted a KEGG pathway enrichment analysis on these target genes of DEMs. As shown in Fig. 6a and b, we identified 18 critical pathways for up-regulated target genes and the top 20 pathways for down-regulated target genes, respectively. KEGG analysis suggested that the up-regulated genes were significantly enriched in insulin resistance, PI3K-Akt signaling pathway, and insulin signaling pathway, while the down-regulated genes mainly participated in the Wnt signaling, cancer, melanogenesis, and proteoglycans in cancer, etc.

Fig. 6
figure 6

KEGG pathway analyses for the target genes of DEMs. a KEGG pathway analysis for up-regulated DEMs. b KEGG pathway analysis for down-regulated DEMs

Construction of DEM-Hub Genes Network

The top 30 hub genes for up-regulated and down-regulated DEMs are shown in Fig. 7a and b, respectively. For up-regulated DEMs, the top 10 hub genes were KRAS, NRAS, CDC42, GDNF, SOS1, PIK3R3, GSK3B, IRS2, GNG12 and PRKCA. For down-regulated DEMs, the top 10 hub genes were NR3C1, PPARGC1A, SUMO1, MEF2C, FOXO3, PPP1CB, MAP2K1, RARA, RHOC, CDC23, and CREB3L2(Table 3). We constructed the candidate miRNA-hub genes regulatory network to precisely explore the molecular mechanism of these DEMs, which was presented in the form of a diagram (Fig. 8).

Fig. 7
figure 7

Identification of the hub genes for DEMs in the PPI network. a PPI network of the top 30 hub genes for up-regulated DEMs. b PPI network of the top 30 hub genes for down-regulated DEMs

Table 3 Top 10 hub genes of the significantly up-regulated and down-regulated DEMs in the PPI network ranked by MCC
Fig. 8
figure 8

Identification of the hub genes for DEMs in the PPI network. Diamond represents genes, ellipse represents miRNAs. Red represents up-regulated DEMs and hub genes, Blue represents down-regulated DEMs and hub genes

Validation of Hub Genes Expression

In this study, 13 hub genes with high interactions such as KRAS, NRAS, CDC42, GDNF, SOS1, PIK3R3, NR3C1, PPARGC1A, SUMO1, MEF2C, FOXO3, PPP1CB and MAP2K1 were verified through GEO dataset GSE114007. Regrettably, GDNF is barely expressed in this database. The expression of KRAS, NRAS, PIK3R3, and SOS1 continue to increase for down-regulated DEMs, which was consistent with the bioinformatics analysis (Fig. 9a-e). For up-regulated DEMs, the expression levels of NR3C1 and PPP1CB were incompatible with the previous screening results; other mRNAs (PPARGC1A, SUMO1, MEF2C, FOXO3 and MAP2K1) were not significantly different(P > 0.05)(Fig. 9f-l). Thus, mir-584-5p-KRAS, mir-183-5p-NRAS, mir-4435-PIK3R3 and mir-4435-SOS1 were identified as four potential regulatory pathways in our study.

Fig. 9
figure 9

The mRNA expression of the top 12 hub genes was determined from the GSE114007 dataset (a). KRAS (b). NRAS (c). CDC42 (d). SOS1 (e). PIK3R3 f). NR3C1 (g). PPARGC1A (h). SUMO1 (i). MEF2C (j). FOXO3 (k). PPP1CB (l). MAP2K1. P < 0.05 was considered statistically significant

Discussion

OA is the most common chronic condition associated with aging and progressive joint dysfunction and the disease with the greatest socioeconomic cost. There is a clear and urgent need to explore the aetiology and progression and support the development of disease-modifying therapies for OA. Our study provided strong evidence for the replication of previous screening integrating RNA-seq with integrated bio-informatics analysis, which may considerably increase the possibility of identifying potential biomarkers and greatly enhance the stability of our findings. To the best of our knowledge, this is the first study to comprehensively explore and identify hub miRNAs of OA from PBL of participants in Southern Chinese population. The study shed a light on plausible etiology of OA, although their biological mechanism still remains to be elucidated in the future.

In the study, we screened out 9 up-regulated DEMs and 10 down-regulated DEMs by exploring overlapped DEMs in the PBL of internal samples and synovial tissue of GSE143514 database. Among those up-regulated DEMs, mir-34a-5p could play an important role in the cell progression of osteosarcoma and could induce chondrocyte apoptosis during OA progression [26]. A previous study reported that mir-138-5p could inhibit the apoptosis and inflammation of chondrocytes and regulate cartilage phenotype and osteogenesis [27]. For down-regulated DEMs, mir-1246 could be involved in the regulation of chondrocyte cell apoptosis and IL-6 secretion from chondrocyte cells [28]. Notably, those critical miRNAs play important roles in regulating biological functions in the process of OA.

Previous studies have shown that TFs can regulate miRNA expression [29]. We found that ERG1 definitely co-exist as the top transcription factor among those up-regulated DEMs and down-regulated DEMs. EGR1 belongs to the EGR family of C2H2-type zinc-finger protein [30]. EGR1 could accelerate cartilage hyper proliferation and degeneration by regulating KLF5 expression and β-catenin signaling pathway [31]. Specificity protein 1(SP1), as a C2H2-type zinc-finger transcription factor, inhibits the proliferation of chondrocytes and promotes the apoptosis of chondrocytes by binding with miR-145 [32]. Eukaryotic Translation Termination Factor 1(E2F1) could regulate the potentials in modulating articular chondrocytes and participate in the regulatory mechanism of knee OA [33]. This suggest that these candidate DEMs play important roles in the pathogenesis of OA.

The expression of 4 genes (KRAS, NRAS, PIK3R3, SOS1) was consistent with that in GSE114007 dataset from the results of 12 hub genes, representing high priority candidates for therapeutic interventions. KRAS, a heterodimer of RAS, leads to chondrocyte growth by stimulating the MAPK/ERK1/2/Cyclin D1 signaling pathway [34]. RAS has been shown to stabilize KRAS levels in the cytoplasmic membrane, aggregating KRAS to amplify the signal output of the mitogen-activated protein kinase (MAPK) pathway [35]. As a member of the pleiotrophin, Midkine (MK) can interact with RAS and further lead to the formation of MK-RAS-KRAS complexes [36]. Then, this action stimulates the downstream MAPK/ErK1/2 signalling pathway, thereby increasing the expression of cyclin D1 and accelerating cell cycle progression [37]. Previous studies reported positive correlations between chondrocyte death by apoptosis and the severity of OA [38, 39]. Moreover, NRAS is an N-RAS oncogene that encodes a membrane protein travelling between the Golgi body and the plasma membrane. Actually, only a small number of studies investigated the role of NRAS in the development of OA. Huang et al. have reported that NRAS is an up-regulated gene in synovial tissue and associated with the MAPK pathway based on pathway analysis [40]. The MAPK pathway is the most important signal transduction system that mediates osteoarthritis cartilage injury, causing a series of reactions, such as the increased expression of matrix MMPs, chondrocyte apoptosis and cartilage destruction [40]. As a key subunit of the IAPI3K protein complex [41], PIK3R3 regulates anabolism effector T cell differentiation, T cell senescence, and activation-induced cell death [42]. PIK3R3 is down-regulated in human cartilage tissue with OA and also mediates the effects of miR-1236 on chondrocyte proliferation and apoptosis [43]. Additionally, SOS1 is a guanine nucleotide exchange factor (GEF) for RAS involved in the RAS signaling pathway [44]. Major components of RAS, including ACE, AT1R, and AT2R, are expressed in synovial tissue in humans and participate in the pathogenesis of OA [45]. Kawakami et al. confirmed that chondrocytes partially expressed AT1R and AT2R and presented increased expression with the stimulation of inflammatory cytokines IL-1 [46].

This study has several merits. First, we have systematically characterized the profile of miRNAs and mRNAs expression between OA patients and healthy controls (blood samples) integrating with public data GSE143514(synovial tissues) to screen overlapped DEMs, which contain different tissues from participants. Second, we performed miRNA-associated integrated network analysis to identify potential hub genes by bioinformatics analysis. More importantly, the acquisition of blood samples is actually feasible as a means of early screening and diagnosis compared with cartilage samples and/or synovial tissues in clinical practice.

Also, there are still some limitations need to be considered. First of all, the sample size of the integrated analysis, neither the internal data nor the GEO dataset, is sufficient to explore and verify those candidate miRNAs and key genes. Considering the veracity and credibility, further studies with larger sample size are still needed to validate the findings. And, although we employed the internal samples and the public databases to comprehensively screen the critical miRNAs and key genes, the involvement of these factors and their value as therapeutic targets need to be validated in studies on their function in vitro and in vivo.

In summary, using integrated analysis, we have detected a series of potential miRNA-target gene pathways (mir-584-5p-KRAS, mir-183-5p-NRAS, mir-4435-PIK3R3, and mir-4435-SOS1) that represent new mediators of abnormal gene expression and high priority targets for therapeutic interventions. Our findings could advance in basic understanding of risk factors and molecular mechanisms of OA disease. Screening from large populations also remains to be further elucidated and validated.

Availability of data and materials

This manuscript contains previously unpublished data. Internal dataset has been uploaded to jianguoyun software (https://www.jianguoyun.com/p/DSPmrX8Q1Iv9CRjW4ZgE). The other datasets analyzed in the present study are available in Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/).

Abbreviations

OA:

Osteoarthritis

GWAS:

Genome-wide association studies

miRNA:

MicroRNA

DEMs:

Differentially expressed microRNAs

PBL:

Peripheral blood leukocytes

GEO:

Gene Expression Omnibus

GO:

Gene Ontology

BP:

Biological process

MF:

Molecular function

CC:

Cellular component

KEGG:

Kyoto Encyclopedia of Genes and Genomes

PPI:

Protein-protein interaction

References

  1. Kulkarni P, Martson A, Vidya R, Chitnavis S, Harsulkar A. Pathophysiological landscape of osteoarthritis. Adv Clin Chem. 2021;100:37–90.

    Article  PubMed  Google Scholar 

  2. Ansari MY, Ahmad N, Haqqi TM. Oxidative stress and inflammation in osteoarthritis pathogenesis: Role of polyphenols. Biomed Pharmacother. 2020;129:110452.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Hawker GA. Osteoarthritis is a serious disease. Clin Exp Rheumatol. 2019;37(Suppl 120(5)):3–6.

    PubMed  Google Scholar 

  4. Zhang Z, Huang C, Jiang Q, Zheng Y, Liu Y, Liu S, et al. Guidelines for the diagnosis and treatment of osteoarthritis in China (2019 edition). Ann Translational Med. 2020;8(19):1213.

    Article  Google Scholar 

  5. Loughlin J. Genetic indicators and susceptibility to osteoarthritis. Br J Sports Med. 2011;45(4):278–82.

    Article  CAS  PubMed  Google Scholar 

  6. Xie F, Liu YL, Chen XY, Li Q, Zhong J, Dai BY, et al. Role of MicroRNA, LncRNA, and Exosomes in the Progression of Osteoarthritis: A Review of Recent Literature. Orthop Surg. 2020;12(3):708–16.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Malemud CJ. MicroRNAs and Osteoarthritis. Cells. 2018;7(8):92.

  8. Cheleschi S, Gallo I, Barbarino M, Giannotti S, Mondanelli N, Giordano A, et al. MicroRNA Mediate Visfatin and Resistin Induction of Oxidative Stress in Human Osteoarthritic Synovial Fibroblasts Via NF-kappaB Pathway. Int J Mol Sci. 2019;20(20):5200.

  9. 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.

    Article  CAS  Google Scholar 

  10. van Meurs JB. Osteoarthritis year in review 2016: genetics, genomics and epigenetics. Osteoarthr Cartil. 2017;25(2):181–9.

    Article  Google Scholar 

  11. Li YH, Tavallaee G, Tokar T, Nakamura A, Sundararajan K, Weston A, et al. Identification of synovial fluid microRNA signature in knee osteoarthritis: differentiating early- and late-stage knee osteoarthritis. Osteoarthr Cartil. 2016;24(9):1577–86.

    Article  CAS  Google Scholar 

  12. Aryal B, Singh AK, Rotllan N, Price N, Fernandez-Hernando C. MicroRNAs and lipid metabolism. Curr Opin Lipidol. 2017;28(3):273–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Tan X, Fu Y, Chen L, Lee W, Lai Y, Rezaei K, et al. miR-671-5p inhibits epithelial-to-mesenchymal transition by downregulating FOXM1 expression in breast cancer. Oncotarget. 2016;7(1):293–307.

    Article  PubMed  Google Scholar 

  14. Nan A, Chen L, Zhang N, Liu Z, Yang T, Wang Z, et al. A novel regulatory network among LncRpa, CircRar1, MiR-671 and apoptotic genes promotes lead-induced neuronal cell apoptosis. Arch Toxicol. 2017;91(4):1671–84.

    Article  CAS  PubMed  Google Scholar 

  15. Nakamura A, Rampersaud YR, Nakamura S, Sharma A, Zeng F, Rossomacha E, et al. microRNA-181a-5p antisense oligonucleotides attenuate osteoarthritis in facet and knee joints. Ann Rheum Dis. 2019;78(1):111–21.

    Article  CAS  PubMed  Google Scholar 

  16. Ntoumou E, Tzetis M, Braoudaki M, Lambrou G, Poulou M, Malizos K, et al. Serum microRNA array analysis identifies miR-140-3p, miR-33b-3p and miR-671-3p as potential osteoarthritis biomarkers involved in metabolic processes. Clin Epigenetics. 2017;9:127.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Chu M, Zhu X, Wang C, Rong J, Wang Y, Wang S, et al. The rs4238326 polymorphism in ALDH1A2 gene potentially associated with non-post traumatic knee osteoarthritis susceptibility: a two-stage population-based study. Osteoarthr Cartil. 2017;25(7):1062–7.

    Article  CAS  Google Scholar 

  18. Langmead B. Aligning short sequencing reads with Bowtie. Curr Protoc Bioinformatics. 2010;Chapter 11:Unit 11 7.

    PubMed  Google Scholar 

  19. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30(7):923–30.

    Article  CAS  PubMed  Google Scholar 

  20. Clough E, Barrett T. The Gene Expression Omnibus Database. Methods Mol Biol. 2016;1418:93–110.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Zhou Y, Wang Z, Chen X, Zhang J, Yang L, Liu S, et al. Identification of differentially expressed miRNAs and mRNAs in synovial of osteoarthritis via RNA-sequencing. BMC Med Genet. 2020;21(1):46.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Fonseka P, Pathan M, Chitti SV, Kang T, Mathivanan S. FunRich enables enrichment analysis of OMICs datasets. J Mol Biol. 2021;433(11):166747.

    Article  CAS  PubMed  Google Scholar 

  23. Kramarz B, Lovering RC. Gene Ontology: A Resource for Analysis and Interpretation of Alzheimer’s Disease Data. In: Wisniewski T, editor. Alzheimer’s Disease. Brisbane (AU): Codon Publications Copyright: The Authors.; 2019.

    Google Scholar 

  24. Du J, Yuan Z, Ma Z, Song J, Xie X, Chen Y. KEGG-PATH: Kyoto encyclopedia of genes and genomes-based pathway analysis using a path analysis model. Mol BioSyst. 2014;10(9):2441–7.

    Article  CAS  PubMed  Google Scholar 

  25. Zhou BF. Predictive values of body mass index and waist circumference for risk factors of certain related diseases in Chinese adults--study on optimal cut-off points of body mass index and waist circumference in Chinese adults. Biomed Environ Sci. 2002;15(1):83–96.

    PubMed  Google Scholar 

  26. Dong C, Wang X, Li N, Zhang K, Wang X, Zhang H, et al. microRNA-mediated GAS1 downregulation promotes the proliferation of synovial fibroblasts by PI3K-Akt signaling in osteoarthritis. Exp Ther Med. 2019;18(6):4273–86.

    CAS  PubMed  PubMed Central  Google Scholar 

  27. Shi J, Cao F, Chang Y, Xin C, Jiang X, Xu J, et al. Long non-coding RNA MCM3AP-AS1 protects chondrocytes ATDC5 and CHON-001 from IL-1β-induced inflammation via regulating miR-138-5p/SIRT1. Bioengineered. 2021;12(1):1445–56.

    Article  CAS  PubMed  Google Scholar 

  28. Qi K, Lin R, Xue C, Liu T, Wang Y, Zhang Y, et al. Long Non-Coding RNA (LncRNA) CAIF is Downregulated in Osteoarthritis and Inhibits LPS-Induced Interleukin 6 (IL-6) Upregulation by Downregulation of MiR-1246. Med Sci Monit. 2019;25:8019–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Pathan M, Keerthikumar S, Ang CS, Gangoda L, Quek CY, Williamson NA, et al. FunRich: An open access standalone functional enrichment and interaction network analysis tool. Proteomics. 2015;15(15):2597–601.

    Article  CAS  PubMed  Google Scholar 

  30. Uddin MN, Li M, Wang X. Identification of Transcriptional Markers and microRNA-mRNA Regulatory Networks in Colon Cancer by Integrative Analysis of mRNA and microRNA Expression Profiles in Colon Tumor Stroma. Cells. 2019;8(9):1054.

  31. Zhang HY, Liu Q, Yang HX, Shi LQ, Wang P, Xie MJ, et al. Early growth response 1 reduction in peripheral blood involving condylar subchondral bone loss. Oral Dis. 2019;25(7):1759–68.

    Article  CAS  PubMed  Google Scholar 

  32. Xue H, Yu P, Wang WZ, Niu YY, Li X. The reduced lncRNA NKILA inhibited proliferation and promoted apoptosis of chondrocytes via miR-145/SP1/NF-κB signaling in human osteoarthritis. Eur Rev Med Pharmacol Sci. 2020;24(2):535–48.

    CAS  PubMed  Google Scholar 

  33. Wang K, Li F, Yuan Y, Shan L, Cui Y, Qu J, et al. Synovial Mesenchymal Stem Cell-Derived EV-Packaged miR-31 Downregulates Histone Demethylase KDM2A to Prevent Knee Osteoarthritis. Mol Ther Nucleic Acids. 2020;22:1078–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Deng Q, Yu X, Deng S, Ye H, Zhang Y, Han W, et al. Midkine promotes articular chondrocyte proliferation through the MK-LRP1-nucleolin signaling pathway. Cell Signal. 2020;65:109423.

    Article  CAS  PubMed  Google Scholar 

  35. Inder KL, Lau C, Loo D, Chaudhary N, Goodall A, Martin S, et al. Nucleophosmin and nucleolin regulate K-Ras plasma membrane interactions and MAPK signal transduction. J Biol Chem. 2009;284(41):28410–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Muramatsu T. Structure and function of midkine as the basis of its pharmacological effects. Br J Pharmacol. 2014;171(4):814–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Cully M, Downward J. SnapShot: Ras Signaling. Cell. 2008;133(7):1292–e1.

    Article  CAS  PubMed  Google Scholar 

  38. Zamli Z, Sharif M. Chondrocyte apoptosis: a cause or consequence of osteoarthritis? Int J Rheum Dis. 2011;14(2):159–66.

    Article  PubMed  Google Scholar 

  39. Mort JS, Billington CJ. Articular cartilage and changes in arthritis: matrix degradation. Arthritis Res. 2001;3(6):337–41.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Huang H, Zheng J, Shen N, Wang G, Zhou G, Fang Y, et al. Identification of pathways and genes associated with synovitis in osteoarthritis using bioinformatics analyses. Sci Rep. 2018;8(1):10050.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Pons S, Asano T, Glasheen E, Miralpeix M, Zhang Y, Fisher TL, et al. The structure and function of p55PIK reveal a new regulatory subunit for phosphatidylinositol 3-kinase. Mol Cell Biol. 1995;15(8):4453–65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Sogkas G, Adriawan IR, Dubrowinskaja N, Atschekzei F, Schmidt RE. Homeostatic and pathogenic roles of PI3Kdelta in the human immune system. Adv Immunol. 2020;146:109–37.

    Article  CAS  PubMed  Google Scholar 

  43. Wang WT, Huang ZP, Sui S, Liu JH, Yu DM, Wang WB. microRNA-1236 promotes chondrocyte apoptosis in osteoarthritis via direct suppression of PIK3R3. Life Sci. 2020;253:117694.

    Article  CAS  PubMed  Google Scholar 

  44. Nimnual A, Bar-Sagi D. The two hats of SOS. Sci STKE. 2002;2002(145):pe36.

    Article  PubMed  Google Scholar 

  45. Hennig A, Markwart R, Wolff K, Schubert K, Cui Y, Prior IA, et al. Feedback activation of neurofibromin terminates growth factor-induced Ras activation. Cell Commun Signal. 2016;14:5.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Kawakami Y, Matsuo K, Murata M, Yudoh K, Nakamura H, Shimizu H, et al. Expression of Angiotensin II Receptor-1 in Human Articular Chondrocytes. Arthritis. 2012;2012:648537.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgments

We acknowledge all participants and staff in the collection of participants.

Funding

This study was supported by Local High-level University (cultivation) project in Shanghai(E1–2602–21-201006-6), and supported by the Shanghai Municipal Health Commission (201840297), and supported by The Training Planned Fund of Academic Leaders from Shanghai Pudong New Area Health System (PWR 12018–09).

Author information

Authors and Affiliations

Authors

Contributions

Ying Jiang drafted the protocol and wrote the final paper. Liying Jiang contributed to interpretation of results and made critical revisions. Yi Shen, Liyan Ding, Shengli Xia participated in the data collection. All authors have reviewed the final version of the manuscript and approved it for publication.

Corresponding author

Correspondence to Liying Jiang.

Ethics declarations

Ethics approval and consent to participate

The study protocol has been approved by the ethics board of Shanghai University of Medicine & Health Sciences. All methods were carried out in accordance with relevant guidelines and informed consent was obtained from all participants.

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

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

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jiang, Y., Shen, Y., Ding, L. et al. Identification of transcription factors and construction of a novel miRNA regulatory network in primary osteoarthritis by integrated analysis. BMC Musculoskelet Disord 22, 1008 (2021). https://doi.org/10.1186/s12891-021-04894-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12891-021-04894-2

Keywords