- Open Access
Identification of transcription factors and construction of a novel miRNA regulatory network in primary osteoarthritis by integrated analysis
BMC Musculoskeletal Disorders volume 22, Article number: 1008 (2021)
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.
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).
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.
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.
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 . 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 . 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 . 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 . Recent breakthroughs have revealed numerous examples of miRNAs involvement in normal development and disease . 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 . And, much attention has recently focused on miRNA because increasing evidence indicates that miRNA affect gene transcription through a number of regulatory processes .
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 . 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 . 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
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 . 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  were used to compare DEMs, and feature Counts  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 . The GEO dataset GSE143514, including five OA patients and three controls were screened . 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 , 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) . 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 . 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.
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), . 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).
A schematic diagram of the study design is displayed in Fig. 1, and the characteristics of the internal database are shown in Table 1.
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.
` 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.
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).
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.
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.
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).
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.
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 . A previous study reported that mir-138-5p could inhibit the apoptosis and inflammation of chondrocytes and regulate cartilage phenotype and osteogenesis . For down-regulated DEMs, mir-1246 could be involved in the regulation of chondrocyte cell apoptosis and IL-6 secretion from chondrocyte cells . 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 . 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 . EGR1 could accelerate cartilage hyper proliferation and degeneration by regulating KLF5 expression and β-catenin signaling pathway . 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 . Eukaryotic Translation Termination Factor 1(E2F1) could regulate the potentials in modulating articular chondrocytes and participate in the regulatory mechanism of knee OA . 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 . 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 . As a member of the pleiotrophin, Midkine (MK) can interact with RAS and further lead to the formation of MK-RAS-KRAS complexes . Then, this action stimulates the downstream MAPK/ErK1/2 signalling pathway, thereby increasing the expression of cyclin D1 and accelerating cell cycle progression . 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 . 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 . As a key subunit of the IAPI3K protein complex , PIK3R3 regulates anabolism effector T cell differentiation, T cell senescence, and activation-induced cell death . PIK3R3 is down-regulated in human cartilage tissue with OA and also mediates the effects of miR-1236 on chondrocyte proliferation and apoptosis . Additionally, SOS1 is a guanine nucleotide exchange factor (GEF) for RAS involved in the RAS signaling pathway . Major components of RAS, including ACE, AT1R, and AT2R, are expressed in synovial tissue in humans and participate in the pathogenesis of OA . Kawakami et al. confirmed that chondrocytes partially expressed AT1R and AT2R and presented increased expression with the stimulation of inflammatory cytokines IL-1 .
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/).
Genome-wide association studies
Differentially expressed microRNAs
Peripheral blood leukocytes
Gene Expression Omnibus
Kyoto Encyclopedia of Genes and Genomes
Kulkarni P, Martson A, Vidya R, Chitnavis S, Harsulkar A. Pathophysiological landscape of osteoarthritis. Adv Clin Chem. 2021;100:37–90.
Ansari MY, Ahmad N, Haqqi TM. Oxidative stress and inflammation in osteoarthritis pathogenesis: Role of polyphenols. Biomed Pharmacother. 2020;129:110452.
Hawker GA. Osteoarthritis is a serious disease. Clin Exp Rheumatol. 2019;37(Suppl 120(5)):3–6.
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.
Loughlin J. Genetic indicators and susceptibility to osteoarthritis. Br J Sports Med. 2011;45(4):278–82.
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.
Malemud CJ. MicroRNAs and Osteoarthritis. Cells. 2018;7(8):92.
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.
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.
van Meurs JB. Osteoarthritis year in review 2016: genetics, genomics and epigenetics. Osteoarthr Cartil. 2017;25(2):181–9.
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.
Aryal B, Singh AK, Rotllan N, Price N, Fernandez-Hernando C. MicroRNAs and lipid metabolism. Curr Opin Lipidol. 2017;28(3):273–80.
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.
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.
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.
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.
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.
Langmead B. Aligning short sequencing reads with Bowtie. Curr Protoc Bioinformatics. 2010;Chapter 11:Unit 11 7.
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.
Clough E, Barrett T. The Gene Expression Omnibus Database. Methods Mol Biol. 2016;1418:93–110.
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.
Fonseka P, Pathan M, Chitti SV, Kang T, Mathivanan S. FunRich enables enrichment analysis of OMICs datasets. J Mol Biol. 2021;433(11):166747.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Muramatsu T. Structure and function of midkine as the basis of its pharmacological effects. Br J Pharmacol. 2014;171(4):814–26.
Cully M, Downward J. SnapShot: Ras Signaling. Cell. 2008;133(7):1292–e1.
Zamli Z, Sharif M. Chondrocyte apoptosis: a cause or consequence of osteoarthritis? Int J Rheum Dis. 2011;14(2):159–66.
Mort JS, Billington CJ. Articular cartilage and changes in arthritis: matrix degradation. Arthritis Res. 2001;3(6):337–41.
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.
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.
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.
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.
Nimnual A, Bar-Sagi D. The two hats of SOS. Sci STKE. 2002;2002(145):pe36.
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.
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.
We acknowledge all participants and staff in the collection of participants.
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).
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
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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.
About this article
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
- Circulating microRNAs
- RNA sequencing
- Bioinformatics analysis