Integrated Analysis of miRNAs and Gene Expression Profiles Reveals Potential Biomarkers for Osteoarthritis

Purpose: Currently, the early diagnosis and treatment of osteoarthritis (OA) remain a challenge. In the present study, we attempted to explore potential biomarkers for the diagnosis and treatment of OA.

Methods: The differentially expressed genes (DEGs) were identified based on three mRNA datasets of synovial tissues for OA patients and normal controls downloaded from the Gene Expression Omnibus (GEO) database. Furthermore, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were used for evaluating gene function related categories. Then, miRNA sequencing was performed for differentially expressed miRNAs’ identification. Finally, weighted gene co-expression network analysis (WGCNA) was performed for genes detected by the three mRNA datasets and a competing endogenous RNA (ceRNA) network with DEGs and differentially expressed microRNAs (miRNAs) was constructed for central genes identification. In addition, the relationship between central gene expression and immune infiltration was analyzed, and the candidate agents for OA were predicted based on the Connectivity Map database. Quantitative RT-PCR (qRT-PCR), Western blotting analysis, and immunofluorescent staining were performed to validate the expression levels of differentially expressed miRNAs and differentially expressed target genes in normal and OA tissues and chondrocytes. MiRNA–mRNA network was also validated in chondrocytes in vitro.

Results: A total of 259 DEGs and 26 differentially expressed miRNAs were identified, among which 94 miRNA–mRNA interactions were predicted. The brown module in WGCNA was most closely correlated with the clinical traits of OA. After overlapping the brown module genes with miRNA–mRNA pairs, 27 miRNA–mRNA pairs were obtained. A ceRNA network was constructed with 5505 lncRNA–miRNA–mRNA interactions. B-cell translocation gene 2(BTG2), Abelson-related gene (ABL2), and vascular endothelial growth factor A (VEGFA) were identified to be the central genes with good predictive performance, which were significantly correlated with immune cell infiltration in OA, reflected by declined activated dendritic cells (aDCs), and elevated contents of B cells, macrophages, neutrophils, and T helper cells. Anisomycin, MG-132, thapsigargin, and lycorine were predicted to be the potential candidate agents for OA intervention. In vitro, the expression levels of differentially expressed miRNAs and biomarkers identified in the present study were consistent with the results obtained in normal or OA knee cartilage tissues and chondrocytes. Furthermore, BTG2 was identified to be negatively regulated by miR-125a-5p.

Conclusion: BTG2, ABL2, and VEGFA can be regarded as potential predictive and treatment biomarkers for OA, which might guide the clinical therapy of OA.

Introduction

Osteoarthritis (OA) is the most common type of arthritis characterized by the loss of articular cartilage under the joints and bones. Its main symptoms include joint pain, swelling, and stiffness, which severely affect physical function and the quality of life. OA typically affects the joints, but it also shows a relationship with cardiovascular diseases, such as myocardial infarction (Schieir et al., 2017) and coronary heart disease (Chung et al., 2016). It is estimated that OA affects 242 million people worldwide and the frequency of OA in women is higher than in men (Breedveld, 2004). OA is the primary cause of joint damage and irreversible disability. Thus, the early diagnosis and aggressive treatment management of OA are necessary for improving the prognosis of this disease.

Emerging evidence has revealed that biomarker identification facilitates the improvement of diagnosis and treatment of OA (Watt, 2018). The level of serum superoxide dismutase 2 (SOD2) is elevated in OA patients indicating OA related oxidative stress and has been suggested as the candidate diagnostic biomarker (Fernandez-Moreno et al., 2011). The baseline collagen type II cleavage neoepitope (C2C) level in serum predicts the risk of cartilage damage in OA (Poole et al., 2016). C-reactive protein (CRP) is one of the inflammatory markers strongly associated with the incidence and progression of OA (Saberi et al., 2016). The circulating lipopolysaccharide associated with systematic inflammation is the indicator for an abundance of activated synovial macrophage abundance and osteophyte severity (Watt, 2018). However, more biomarkers are indeed needed to more accurately predict the pathological process of OA so as to achieve the goal of multi-targeted therapy for OA.

The development of next-generation technologies such as large-scale sequencing and multianalyte assay has recently facilitated the discovery of a slew of potential biomarkers for OA (Cong et al., 2017; Xu et al., 2019; Zhang et al., 2019; Hong et al., 2020). For example, Li et al. (2017), identified several differentially expressed genes, including a disintegrin and metalloproteinase 8 (Adam8), arginase 1 (Arg1) and microRNAs (miRNAs) such as miRNA-200b, which might play key roles in the pathogenesis of mice OA based on microarray analysis. Furthermore, functional disruption of collagen in cartilage has been shown to be a very important contributing factor to OA. Quantitative detection of collagen synthesis may facilitate early identification and prediction of OA (Zhang et al., 2019). However, none of these biomarkers is specific, which makes diagnosis and treatment of OA still a challenge. In the present study, we added our own miRNA microarray data and combined them with gene expression datasets deposited in the Gene Expression Omnibus (GEO) database to perform a comprehensive bioinformatic analysis of potential biomarkers for OA. Furthermore, a new ceRNA network was constructed, and drug prediction analysis for OA was performed. In addition, select differentially expressed miRNAs and their target genes were further validated by in vitro experiments. The promising potential biomarkers with predictive function in the pathogenesis and treatment of OA identified in this study are expected to guide novel interventions of OA.

Methods

Data Downloads

Three mRNA datasets of joint synovial samples [accession numbers GSE 1919 (Ungethuem et al., 2010), GSE55235 (Woetzel et al., 2014), and GSE55457 (Woetzel et al., 2014) ] were available from NCBI GEO database. The dataset of GSE1919 included joint synovial tissues of five OA patients and five normal donors. The downloaded GSE55235 dataset contained joint synovial tissue samples from ten OA patients and ten healthy joints of individuals who suffered from fatal accidents. From the GSE55457 dataset, the mRNA data of synovial membranes from ten OA patients and ten normal donors were retrieved. The txt format series matrix files were downloaded, and the batch effect was removed by the surrogate variable analysis (SVA) algorithm. The differentially expressed genes (DEGs) between OA and normal controls were identified by the limma package (version 3.46.0) (Smyth et al., 2005; Ritchie et al., 2015). p < 0.05 and |logFC(fold change)|>1 were set as the cutoff value. The key steps involved in the bioinformatic analysis of datasets obtained from GEO are illustrated in Figure 1.

FIGURE 1

FIGURE 1. A flow chart showing the steps involved in the bioinformatic analysis of datasets obtained from GEO. The differentially expressed genes (DEGs) were identified based on three mRNA datasets. miRNA sequencing was performed for the identification of differentially expressed miRNAs. Weighted gene co-expression network analysis (WGCNA) was performed for genes detected using three mRNA datasets, and a ceRNA network with DEGs and differentially expressed microRNAs (miRNAs) was constructed for central genes identification. The relationship between central gene expression and immune infiltration was analyzed, and the candidate agents or small molecules for potential OA therapy were predicted based on the Connectivity Map database.

Function Enrichment Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were used for evaluating gene function related categories. DEGs were annotated based on GO and KEGG databases using the ClusterProfiler package in R (version 3.18.1) (Yu et al., 2012). The significantly enriched GO and pathway terms were screened with p < 0.05 and q < 0.05.

MiRNA Sequencing and Differential Expression Analysis

Joint synovial tissues from OA patients (five pairs) and normal controls (five pairs) were obtained under surgery in Guangdong Provincial Hospital of Traditional Chinese Medicine (NO. BF 2021-004). Approval was obtained from the ethics committee of Guangdong Provincial Hospital of Traditional Chinese Medicine, and the study procedures complied with the declaration of Helsinki (World Medical Association, 2013). Total RNA was prepared from synovial tissues and purified using the mirVana miRNA isolation kit (Ambion, Austin, TX, United States) according to the manufacturer’s protocols. The quantified RNA was ligated with 3′ and 5′ adaptors. The miRNA-adaptor ligation products were subjected to qRT-PCR amplification. Then, the cDNA library was sequenced on the DNBSEQ platform (Beijing Genomic Institute, Shenzhen, China). The clean reads were aligned to the human reference genome (the Human Genome U95 Set and Human Genome U133A), and the differentially expressed miRNAs were investigated based on the FPKM method by using the limma package. A p < 0.05 was considered as statistically significant.

Weighted Gene Co-Expression Network Analysis

WGCNA facilitates performing weighted correlation network analysis and describing the correlation structure between gene expression profiles, image data, genetic marker data, and proteomics data (Langfelder and Horvath, 2008). The gene co-expression network for DEGs in GSE 1919, GSE55235, and GSE55457 datasets were generated by the WGCNA-R package. The soft threshold was determined using the “sft$powerEstimate” function, and the topological overlap matrix (TOM) was obtained from the adjacency matrix. Hierarchical clustering was performed for TOM to evaluate the structure of the cluster tree. Different branches (colors) of the cluster tree represented different modules.

Competitive Endogenous RNA Network Construction

The differentially expressed miRNA interacted genes were predicted by miRDB (http://mirdb.org/) (Chen and Wang, 2020), miRTarBase (https://mirtarbase.cuhk.edu.cn/∼miRTar Base/miRTarBase_2022/php/index.php) (Huang et al., 2022), and TargetScan (https://www.targetscan.org/vert_72/) (Agarwal et al., 2015) databases. The DEGs commonly predicted by the three databases were selected for further analysis. The long non-coding RNA (lncRNA)-miRNA interactions were predicted by DIANA-LncBase database (www.microrna.gr/LncBase) (Karagkouni et al., 2020). Based on the correlations between lncRNAs, miRNAs, and mRNAs, a ceRNA network was constructed using Cytoscape software (version3.8.2) (Kohl et al., 2011; Ma et al., 2018). Critical genes were analyzed by connectivity degrees, and the diagnostic performance of interesting genes was validated by receiver operating characteristic (ROC) curve analysis.

Correlation Analysis of Critical Genes With Immune Infiltration

Based on gene expression profiles of different OA samples, the proportions of 28 types of immune cells were deduced by the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm and the relative contents of immune cells were depicted by the “vioplot” package (version 0.3.7) (Hu, 2020). The interplay between different types of immune cells was analyzed by the “corrplot” package (version 0.92) (Zhang et al., 2021). The correlation between gene expression level and immune cell proportion was investigated using Spearman correlation analysis. A p < 0.05 was considered as statistically significant.

Drug Small Molecule Prediction

The Connectivity Map (CMap) is a publicly available database that collects gene expression profiles perturbed by a series of small molecules in various cell lines (Langfelder and Horvath, 2008). DEGs were divided into upregulated and downregulated ones. The expression profiles of upregulated and downregulated DEGs were compared with those in the CMap database. The relative correlations between query signature and reference gene expression profiles of individual drugs in the CMap database are output as enrichment scores ranging from −1 to 1. When the enrichment scores are closed to 1, it means that the corresponding small molecules can contribute to the progression of the disease; otherwise, the small molecules can reverse the gene signals induced by the disease. Thus, the small molecules with enrichment scores closed to −1 and p < 0.05 were selected as the potential candidate drugs for OA therapy.

Gene Set Enrichment Analysis

GSEA was developed as a method for gene expression profile analysis (Subramanian et al., 2005). The gene sets were pre-ranked based on the differential expression between OA and normal control groups. The ordered gene sets were subjected to GSEA. The number of permutations was set as 1000, and the permutation type was set as phenotype.

Chondrocyte Isolation and Culture

This study was performed in accordance with the animal ethics standards and with the approval of the animal ethics committee at the Guangzhou University of Chinese Medicine. Chondrocytes were isolated from embryonic 2-week-old male Sprague–Dawley (SD) rats (Laboratory Animal Center of Guangzhou University of Chinese Medicine, Guangzhou, China) according to an established protocol (Zhang et al., 2014). In brief, cartilage tissue samples were trimmed into small pieces and digested in pronase in PBS (2 mg/ml) for an hour and then in DMEM medium (Gibco, Life Technologies, Carlsbad, CA, United States) containing collagenase D (3 mg/ml) (Roche) or collagenase type Ⅱ (2 mg/ml) (Gibco) for 3 h at 37°C in an incubation shaker. Dissociated cells were cultured in DMEM/F12 medium containing 10% fetal bovine serum (FBS; Gibco, Life Technologies, Carlsbad, CA, United States), 100 U ml⁻¹ of penicillin, and 100 mg ml⁻¹ of streptomycin (Gibco, Life Technologies, Carlsbad, CA, United States). The cells were cultured in DMEM/F12 at 37°C, 5% CO_2, and saturated humidity. The medium was changed every 2 days, and non-adherent cells were discarded. All chondrocytes used in this study were between passages 2 and 4. To establish an OA chondrocyte model, chondrocytes were treated with DMEM/F12 containing interleukin (IL)-1β (10 ng/ml, Sigma-Aldrich, St Louis, Missouri, United States).

Quantitative Real-Time Polymerase Chain Reaction

The qRT-PCR was performed to evaluate the expression levels of differentially expressed miRNAs, including miR-424-5p, miR-125a-5p, and miR-125b-5p and differentially expressed target genes, including BTG2, ABL2, and VEGFA in chondrocytes and knee cartilage tissues of normal and OA patients. OA patients (n = 15) and normal individuals (n = 15) were included in this study. The tissues from subjects were obtained under surgery with informed consent. This study was approved by the Ethics Committee of Guangdong Provincial Hospital of Traditional Chinese Medicine (NO.BF 2021-004).

The total RNA was extracted from chondrocytes or knee cartilage tissues using an RNAiso Plus Kit (Takara Biotechnology, Dalian, Liaoning, China). The PrimeScript reverse transcription reagent kit (Takara Biotechnology, Dalian, Liaoning, China) was applied for reverse transcription. The qRT-PCR was conducted using a SteponePlus (Applied Biosystems, Foster, CA, United States) PCR instrument. SYBR^® Green Premix Pro Taq HS qPCR kit Ⅱ (Accurate Biotechnology, Guangzhou, Guangdong, China) was used for qPCR analysis. The PCR reaction volumes contained 2 μl of cDNA solution, 10 μl of 2 × SYBR^® Green Pro Taq HS Premix Ⅱ, 0.8 μl of PCR forward primer and 0.8 μl of reverse primer (both at 10 μM) and RNase free water up to 20 μl. With U6 selected as an internal reference, the mRNA expression was detected using the 2^−ΔΔCt method. The primer sequences are shown in Table 1. The universal primers in the kit are used as miRNA downstream primers.

TABLE 1

TABLE 1. Primer sequence list.

Western Blotting

The Western blot analysis was performed to detect the protein expression of differentially expressed target genes, including BTG2, ABL2, and VEGFA, in chondrocytes and knee cartilage tissues of normal and OA patients. Total protein from each sample was extracted using Radio Immunoprecipitation Assay (RIPA) lysis buffer (Gibco, Grand Island, NY, United States). Subsequently, 20 μg total proteins from each sample was resolved by 10% SDS-PAGE and transferred onto polyvinylidene difluoride (PVDF, Millipore, Billerica, MA, United States) membranes. After being blocked with 5% skim milk for 90 min at room temperature (RT), the membranes were then incubated with the primary antibody solution overnight at 4°C. Subsequently, the membranes were incubated with a horseradish peroxidase (HRP)–conjugated secondary antibody (1:1000) for 90 min at RT. The protein bands were visualized via a ChemiDoc™ MP Imaging System (Bio-Rad). GAPDH was used as an internal reference. Primary antibodies used in this experiment were from the following source: Rabbit anti-BTG2 (1:1000), mouse anti-ABL2 (1:1000), rabbit anti-VEGFA (1:1000), and GAPDH (1:5000). All primary antibodies were purchased from Invitrogen (Invitrogen inc., Carlsbad, CA, United States).

Immunofluorescent Staining

For immunofluorescence, chondrocytes were fixed with 4% PFA for 30 min and permeabilized in 0.3% Triton X-100 at RT for 30 min. Subsequently, the cells were blocked with 10% normal horse serum at RT for 1 h and then incubated with primary antibodies overnight at 4°C. Next, the primary antibodies were conjugated with Alexa Fluor 488/568 goat anti-rabbit or mouse antibodies (1:300) for identification, and then the nucleus was counter-stained with DAPI. Finally, the slides were mounted and observed under a fluorescence microscope (Olympus IX73). Primary antibodies used in this study were as follows: Rabbit anti-BTG2 (1:1000), mouse anti-ABL2 (1:1000), and rabbit anti-VEGFA (1:1000). All antibodies were purchased from Invitrogen (Invitrogen inc., Carlsbad, CA, United States).

Luciferase Reporter Assay and Chondrocytes Transfection

The target gene of miR-125a-5p was predicted by the online software Targetscan 7.1 (www.targetscan.org/vert_71/). DNA fragments containing the miR-125a-5p binding site in the 3′-untranslated regions (3′-UTR) of BTG2 and fragments containing the mutant miR-125a-5p binding site were separately inserted into the luciferase reporter gene plasmid to obtain the BTG2-WT (wild type) and BTG2-MUT (mutant type) plasmids, respectively. Next, the miR-125a-5p mimics and BTG2-WT/BTG2-MUT plasmids were co-transfected into 293T cells, and the luciferase activities of the transfected cells were detected using luciferase assay. In addition, 293T cells were co-transfected with negative control (NC) oligonucleotide and BTG2-WT/BTG2-MUT plasmids as control. Chondrocytes were seeded onto a six-well plate at a density of 2 × 10⁵ cells/well 2 days prior to transfection. When growing at a confluence of 70%–80%, chondrocytes were transfected with miR-125a-5p mimics, inhibitor, and their negative control (RIBO Biotech Company, Guangzhou, Guangdong, China) using the Lipofectamine 2000 reagent, respectively. At 48 h post-transfection, Western blotting and qRT-PCR were performed to evaluate the effect of cell transfection.

Statistical Analysis

The SPSS 20.0 software (SPSS, Chicago, IL, United States) was applied for statistical analyses. Results are presented as means ± standard deviations (SDs). Multiple groups were compared using one-way ANOVA, while the differences between two groups were analyzed using an unpaired Student’s t-test. A p < 0.05 was considered statistically significant.

Results

Identification of Differentially Expressed Genes of Synovial Tissues of Osteoarthritis Patients Compared With Normal Controls

Three gene expression datasets of synovial tissues were downloaded from the GEO database. For data integration, the batch effects were corrected by surrogate variable analysis (SVA). Before batch effects removal, the samples in different batches were separately clearly, while SVA adjustment brought the samples to a single batch (Figure 2A), indicating the batch effects were effectively removed. Then, the DEGs of three groups of OA patients compared with normal controls were identified by the limma package. Results showed that 259 genes were differentially expressed, including 141 upregulated genes and 118 downregulated genes (Figure 2B). Functional analysis indicated that DEGs were closely related to leukocyte chemotaxis, proliferation, migration, and adhesion related biological processes and significantly enriched in tumor necrosis factor alpha (TNF) signaling and Interleukin (IL)-17 signaling pathways (Figures 2C and D).

FIGURE 2

FIGURE 2. Identification of DEGs in synovial tissues of OA patients compared with normal controls. (A) The samples distribution before and after batch effects removal. PCA, principal component analysis. (B) Volcano plot for 259 DEGs based on three groups of mRNA datasets. (C) GO enrichment analysis of DEGs. (D) The significantly enriched pathways by DEGs.

Competing Endogenous RNA Network and Identification of Critical Genes

By miRNA sequencing, we obtained 26 differentially expressed miRNAs (Figure 3A). According to TargetScan, miRDB, and miRTarBase, there were 1219 miRNA–mRNA interactions mapped to differentially expressed miRNAs. After overlapped with DEGs, 94 miRNA–mRNA interaction pairs were collected. The genes detected in three datasets were subjected to WGCNA. The weight co-expression network was classified into 16 modules, including black (n = 194), blue (n = 957), brown (n = 923), cyan (n = 182), green (n = 222), green-yellow (n = 131), grey (n = 4), grey60 (n = 65), light cyan (n = 80), magenta (n = 273), pink (n = 181), purple (n = 155), red (n = 198), tan (n = 126), turquoise (n = 1021), and yellow (n = 288). The brown module was most significantly correlated with clinical traits (cor = 0.83, p = 2e−13) (Figure 3B). Thus, the genes in the brown module were overlapped with 94 miRNA–mRNA pairs and finally 27 miRNA–mRNA interactions were obtained. In parallel, 452,052 lncRNA–miRNA interactions were predicted by DIANA-LncBase database. Finally, a ceRNA network was constructed with 5,505 lncRNA–miRNA–mRNA interactions (Figure 3C). B-cell translocation gene 2 (BTG2), Abelson-related gene (ABL2), and vascular endothelial growth factor A (VEGFA) were the key mRNAs with highest node degrees and showed relatively good predictive value for OA [BTG2-AUC: 0.902 (0.811–0.994), ABL2-AUC: 0.872 (0.769–0.975), and VEGFA-AUC: 0.973 (0.922–1.000)] (Figures 3D–F).

FIGURE 3

FIGURE 3. ceRNA network construction and the central genes. (A) Heatmap of the differentially expressed miRNAs. The expression profiles of miRNAs were significantly different between OA and control samples. (B) The relationships between WGCNA modules and clinical traits. (C) ceRNA network constructed using 10 differentially expressed miRNAs and 19 DEGs. Red circles represent DEGs, and green squares represent differentially expressed miRNAs. (D–F) The ROC curves for the predictive performance of ABL2, BTG2, and VEGFA for OA.

The Correlation Between Critical Gene Expression and Immune Infiltration

The immune cell content in each sample was predicted (Figure 4A), and the correlation between different types of immune cells was depicted in Figure 4B. In addition, the content of activated dendritic cells (aDCs) was remarkably declined in OA patients, while the contents of B cells, check-point, cytolytic activity, HLA, inflammation-promoting, macrophages, mast cells, neutrophils, parainflammation, T cell co-inhibition, T cell co-stimulation, T helper cells, tumor-infiltrating lymphocytes (TIL), Type I IFN response, and Type II IFN response were obviously elevated (Figure 4C). The expression of BTG2, ABL2, and VEGFA was mostly correlated with the different types of immune cells (Figure 4D).

FIGURE 4

FIGURE 4. The relationship between critical gene expression and immune infiltration. (A) The heatmap of the content for different types of immune cells in OA and control samples. (B) The relationship between different immune cells. (C) The difference in immune cell content between OA and normal samples. (D) The correlation of VEGFA, BTG2, and ABL2 with the content of different immune cells.

Candidate Small Molecules Identification

The DEGs induced expression profiles in OA were compared with that perturbed by small molecules in the CMap database. The top 20 small molecules that showed a significant correlation with query signatures are listed in Table 2. Anisomycin, MG-132, thapsigargin, and lycorine induced expression profiles elicited a highly negative correlation with that of DEGs (Figure 5), suggesting that these molecules may be considered as potential candidates for the intervention of OA and as yet deserve further investigations.

TABLE 2

TABLE 2. Top 20 small molecules predicted by CMap database.

FIGURE 5

FIGURE 5. Molecular structure of the candidate agents for OA. (A–D) The molecular structure for anisomycin, MG-132, thapsigargin, and lycorine, respectively. The DEGs induced expression profiles in OA were compared with that perturbed by small molecules in the CMap database. Anisomycin (A), MG-132 (B), thapsigargin (C), and lycorine (D) induced expression profiles elicited a highly negative correlation with that of the DEGs.

Significant Pathways Involved With Critical Genes

We next examined the major signaling pathways in association with the critical genes identified earlier. Based on GSEA, BTG2 was mainly involved in KEGG spliceosome and adipocytokine signaling pathways. ABL2 was significantly enriched in KEGG P53 signaling and prion diseases related pathways. VEGFA was closely related to KEGG apoptosis and spliceosome pathways (Figure 6).

FIGURE 6

FIGURE 6. Significant pathways implicated with BTG2, ABL2, and VEGFA. Based on GSEA, BTG2 was mainly involved in KEGG spliceosome and adipocytokine signaling pathways. ABL2 was significantly enriched in KEGG P53 signaling and prion diseases related pathways. VEGFA was closely related to KEGG apoptosis and spliceosome pathways.

Expression Levels of Differentially Expressed miRNAs and Target Genes in Normal and Osteoarthritis Tissues and Chondrocytes

QRT-PCR, Western blotting analysis, and immunofluorescent staining were performed to validate the expression levels of differentially expressed miRNAs, including miR-424-5p, miR-125a-5p, and miR-125b-5p and differentially expressed target genes, including BTG2, ABL2, and VEGFA in the ceRNA network. The qRT-PCR results showed that the expression levels of miR-424-5p were significantly downregulated, whereas the expression levels of miR-125a-5p and miR-125b-5p were remarkably upregulated in OA tissues compared with the normal group (Figure 7A). These results were consistent with the miRNA sequencing data. Furthermore, the mRNA expression levels of the target genes, including BTG2, ABL2, and VEGFA in OA tissues were all significantly decreased compared with the normal control group (Figure 7B). Similar effects on the expression of target genes were observed at protein levels. Western blotting analysis revealed an obviously decreased expression of BTG2, ABL2, and VEGFA in BMSCs in OA tissues (Figure 7C).

FIGURE 7

FIGURE 7. The validation of differentially expressed miRNAs and potential biomarkers in normal and OA knee cartilage tissues. (A) QRT-PCR was performed to detect the mRNA expression of miR-424-5p, miR-125a-5p, and miR-125b-5p in normal and OA knee cartilage tissues. (B) QRT-PCR was performed to detect the mRNA expression of BTG2, ABL2, and VEGFA in normal and OA knee cartilage tissues. (C) Western blotting analysis and relative quantification showing the protein expression of BTG2, ABL2, and VEGFA in normal and OA knee cartilage tissues. Values are mean ± SD, *p < 0.05 compared with the Normal 1 group; #p < 0.05 compared with the Normal 2 group.

The immunofluorescent staining results showed that the BTG2, ABL2, and VEGFA positive cell numbers were significantly decreased in the IL-1β induced OA chondrocytes in comparison to that of the control chondrocytes (Figure 8). Taken together, these results indicated that the abovementioned biomarkers exhibit consistent expression trends in OA conditions compared with that of the normal controls.

FIGURE 8

FIGURE 8. The validation of differentially expressed potential biomarkers in normal chondrocytes and IL-1β induced OA chondrocytes. Bright fields and immunofluorescent staining show the protein expression of potential biomarkers, including BTG2, ABL2, and VEGFA positive cells in normal chondrocytes and IL-1β induced OA chondrocytes. The percentage of BTG2, ABL2, and VEGFA positive cells of chondrocytes were also counted and evaluated. Values are mean ± SD. *p < 0.05, compared with the Control 1 group; #p < 0.05 compared with the Control 2 group.

MiR-125a-5p Negatively Regulates B-Cell Translocation Gene 2

The TargetScan software identified that miR-125a-5p could bind with the target sites in the 3′-UTR of BTG2, suggesting that BTG2 might be the target gene of miR-125a-5p (Figure 9A). To determine whether BTG2 was transcriptionally regulated by miR-125a-5p, the WT and MUT BTG2 3′-UTR segment containing the potential target sites for miR-125a-5p was constructed and cloned into the downstream region of a dual-luciferase reporter gene. MiR-125a-5p mimics had no effect on the luciferase activity of 293T cells with BTG2-MUT. However, the luciferase activity of the cells with NC + BTG2-WT was significantly higher than that of the cells with miR-125a-5p mimics + BTG2-WT. Thus, the miR-125a-5p mimics decreased the luciferase activity of cells with BTG2-WT (Figure 9B). These results suggested that miR-125a-5p could directly target BTG2.

FIGURE 9

FIGURE 9. MiR-125a-5p negatively regulates BTG2. (A) The target gene BTG2 of miR-125a-5p was identified using the Targetscan software. (B) The relative luciferase activity between cells was compared after co-transfection of BTG2-WT/BTG2-MUT vectors with miR-125a-5p mimics and negative control. (C) QRT-PCR was performed to detect the mRNA expression of BTG2 in chondrocytes transfected with miR-125a-5p mimics, inhibitor, and their negative control. (D and E) Western blotting analysis and relative quantification showing the protein expression of BTG2 in chondrocytes transfected with miR-125a-5p mimics, inhibitor, and their negative control. GAPDH was used as the internal control. Values are mean ± SD. *p < 0.05.

To further elucidate the effect of miR-125a-5p on BTG2, chondrocytes were transfected with miR-125a-5p mimics, inhibitor, and their negative controls, respectively. The qRT-PCR analysis revealed a significant reduction of BTG2 mRNA level in chondrocytes treated with miR-125a-5p mimics compared with the control; conversely, miR-125a-5p inhibitor obviously enhanced BTG2 mRNA expression (Figure 9C). Similar effects on the expression of BTG2 were observed at protein levels as detected by Western blotting analysis (Figures 9D and E). These results suggested that BTG2 is negatively regulated by miR-125a-5p.

Discussion

Taken together, our results demonstrate that BTG2, ABL2, and VEGFA were identified to be the central genes with good predictive performance, which were significantly correlated with immune cell infiltration in OA, reflected by declined aDCs and elevated contents of B cells, macrophages, neutrophils, and T helper cells. Moreover, anisomycin, MG-132, thapsigargin, and lycorine were predicted to be the candidate agents for potential treatments of OA awaiting further investigations. In vitro, the expression levels of differentially expressed miRNAs and potential biomarkers identified in the present study were consistent with the results obtained in normal or OA knee cartilage tissues and in chondrocytes. In addition, miR-125a-5p negatively regulated BTG2, suggesting the validation of miRNA–mRNA networks.

OA is a chronic degenerative joint disease which is frequently associated with severe joint damage and irreversible functional impairments. OA has affected 52.5 million Americans and generates 185 million economic burden each year (Metsios et al., 2008; Eakin et al., 2017; Pundole and Suarez-Almazor, 2020). It is estimated that there are around 654⋅1 million individuals (40 years and older) with knee OA in 2020 worldwide (Cui et al., 2020). OA has posed a threat to the quality of life of people and has emerged as a public health concern globally. However, at present, the early diagnosis and treatment of OA remain to be an important challenge. Various biomarkers show benefits in improving diagnosis, monitoring disease progression, and predicting prognosis. For the potential of guiding clinical therapies, biomarker identification for OA has attracted increasing attention. Over the last decade, the rise of epigenetics contributed to revealing the pathogenesis of autoimmune disorders and identifying promising candidate biomarkers (Shen et al., 2017; Oliviero et al., 2019; Tavallaee et al., 2020; Wan et al., 2021). Recent evidence has shown that miR-146a-5p is highly expressed in cartilage tissues and serum of OA patients when compared with that of the normal controls and is suggested as a biomarker for OA (Skrzypa et al., 2019). LncRNA PCGEM1 prostate-specific transcript (PCGEM1) is remarkably upregulated in late-stage OA compared with early OA, which has been suggested as an indicator for stratifying OA (Zhao and Xu, 2018). A recent bioinformatics analysis has revealed that the innate immunity-related gene of toll-like receptor 7 (TLR7) is differentially expressed in the synovial membrane, blood, and articular cartilage of OA patients compared with that of the normal controls, suggesting that TLR7 may serve as a potential biomarker for the early diagnosis and treatment of OA (Wang et al., 2020).

CeRNA network containing lncRNA–miRNA–mRNA interactions has been theorized to play an indispensable role in the progression of diseases and could help in exploring the pathogenesis of diseases and biomarker discovery. However, the study of the ceRNA network associated with OA is lacking. In this study, we constructed a ceRNA network and identified BTG2, ABL2, and VEGFA as the critical genes. BTG2, ABL2, and VEGFA were identified to be differentially expressed in synovial tissues of OA patients compared with normal controls based on three groups of mRNA datasets. WGCNA provides information on critical gene function modules and biomarkers for diseases (Tian et al., 2020). The above three genes were in the brown module that was significantly correlated with the clinical traits of OA. All these suggested that BTG2, ABL2, and VEGFA might play a key role in the pathogenesis of OA. In order to validate the critical role of BTG2, ABL2, and VEGFA in OA, we performed miRNA sequencing of synovial tissues of OA patients and normal controls. A total of 26 miRNAs were found to be differentially expressed in OA and the expression profile of miRNAs could clearly distinguish OA samples and normal control samples, indicating that the differentially expressed miRNAs might play important roles in OA. Then, with the mean of a bioinformatic tool for predicting miRNA–mRNA and lncRNA–miRNA interactions, a ceRNA network with DEGs in the brown module and miRNAs with differential expression was constructed. The central role of BTG2, ABL2, and VEGFA was highlighted based on the high node degree.

VEGFA, named vascular endothelial growth factor A, has been reported to be in association with OA (Guan et al., 2020). Zupan et al. (2018) reported that the high expression of VEGFA indicating cartilage or bone degeneration is associated with pre-osteoarthritis. Furthermore, overexpression of miR-140-5p alleviates knee osteoarthritis by downregulating VEGFA in a rat model (Liu et al., 2022). BTG2 was found to be a critical gene related to inflammation in OA. Recent evidence showed that BTG2 is a hub node in the ingenuity pathway network of facet joint OA, which proposed its function in mediating signaling pathways by affecting other genes (Chen et al., 2020). In addition, our ROC analysis indicated that BTG2, ABL2, and VEGFA showed a good potential predictive effect on OA. Thus, the three genes may serve as promising potential predictive biomarkers for OA.

It is well known that synovial inflammation plays a key role in the initiation of OA (de Lange-Brokaar et al., 2012). Recent evidence has shown that there is increased M1 macrophage infiltration and declined mast cell and neutrophil infiltration in OA samples (Yuan et al., 2021). The cartilage damage and synovitis in OA are characterized by abnormal cell proliferation and immune cell infiltration, including T cell and activated macrophages (Lopes et al., 2017). In this study, we found that the aDCs were significantly declined, while B cells, macrophages, T cell co-inhibition, T cell co-stimulation, and T helper cells were significantly elevated in the microenvironment of synovial tissues, which was agreed with the previous findings (Penkava et al., 2020; Van Raemdonck et al., 2020; Zhao et al., 2021). In addition, the expressions of BTG2, ABL2, and VEGFA were significantly correlated with the immune cell infiltration. Thus, the three critical genes may play key roles in the progression of OA by modulating the microenvironment of the articular cartilage.

Furthermore, thapsigargin was predicted to be a candidate small molecule that could reverse the OA process by targeting DEGs. It is reported that insufficient apoptosis of inflammatory cells contributes to the progression of OA, and enhanced apoptosis of chronic inflammatory cells in joints might be the therapeutic target for OA (Bar-Yehuda et al., 2009). Thapsigargin, as a well-known endoplasmic reticular Ca2+ ATPase inhibitor, has been widely used as an apoptosis inducer in tumor cells (Muthukkumar et al., 1995; Chen et al., 2016). Thus, thapsigargin could be suggested as the candidate drug for OA to enhance apoptosis in joints by targeting DEGs, including BTG2, ABL2, and VEGFA. Taken together, the expression levels of differentially expressed miRNAs, including miR-424-5p, miR-125a-5p, and miR-125b-5p and potential biomarkers, including BTG2, ABL2, and VEGFA were validated in normal or OA knee cartilage tissues and chondrocytes. The experimentally validated results are consistent with the bioinformatics analysis. Furthermore, in vitro experiments showed that miR-125a-5p negatively regulated BTG2, suggesting the existence of miRNA–mRNA networks in the pathogenesis of OA.

Limitations

There were several limitations of our study. First, there is no in vitro or in vivo evidence in this study that the screened drugs are useful for the treatment of OA, but there is some literature evidence. Second, validation of potential biomarkers identified in the present study was not performed in vivo models of OA. Further experiments to validate the identified potential biomarkers and illustrate the pathological process and therapeutic targets of OA are therefore necessary.

Conclusion

In summary, BTG2, ABL2, and VEGFA were differentially expressed genes in tissues of OA patients. BTG2, ABL2, and VEGFA were the central genes in the ceRNA network and played key roles in modulating immune cell infiltration in synovial tissues of OA. Moreover, anisomycin, MG-132, thapsigargin, and lycorine were predicted to be the candidate agents for OA. In vitro, the expression levels of differentially expressed miRNAs and potential biomarkers identified in the present study were consistent with the results validated in normal or OA knee tissues and in chondrocytes. BTG2, ABL2, and VEGFA could be regarded as potential diagnostic biomarkers for OA, which might guide the clinical therapy of OA.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding authors.

Ethics Statement

The studies involving human participants were reviewed and approved by the Guangdong Provincial Hospital of Traditional Chinese Medicine (NO.BF 2021-004). The patients/participants provided their written informed consent to participate in this study.

Author Contributions

XC and CW conceived and designed this study. ZHL and ZC wrote the manuscript. ZHL, ZC, and XW acquired and analyzed the main data. HS, ZEL, JW, and XZ collected the samples and contributed analysis tools. All authors read and approved the final manuscript.

Funding

This study was supported by a Translational Collaborative Innovation Project (NO. YN2018HK04) from the Guangdong Provincial Hospital of Traditional Chinese Medicine and the School of Biomedical Sciences, The Chinese University of Hong Kong.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.814645/full#supplementary-material

Abbreviations

OA, osteoarthritis; SOD2, serum superoxide dismutase 2; C2C, collagen type II cleavage neoepitope; CRP, C-reactive protein; Adam8, a disintegrin and metalloproteinase 8; Arg1, arginase 1; miRNAs, microRNAs; GEO, Gene Expression Omnibus; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; WGCNA, weighted gene co-expression network analysis; TOM, topological overlap matrix; lncRNA, long non-coding RNA; ROC, receiver operating characteristic; CMap, connectivity map; GSEA, gene set enrichment analysis; SVA, surrogate variable analysis; TNF, tumor necrosis factor alpha; IL-17, Interleukin-17; ceRNA, competing endogenous RNA; BTG2, B-cell translocation gene 2; ABL2, Abelson-related gene; VEGFA, vascular endothelial growth factor A; aDCs, activated dendritic cells; PCGEM1, PCGEM1 prostate-specific transcript; TLR7, toll-like receptor 7.

References

Agarwal, V., Bell, G. W., Nam, J. W., and Bartel, D. P. (2015). Predicting Effective microRNA Target Sites in Mammalian mRNAs. Elife 4, e05005. doi:10.7554/elife.05005

PubMed Abstract | CrossRef Full Text | Google Scholar

Bar-Yehuda, S., Rath-Wolfson, L., Del Valle, L., Ochaion, A., Cohen, S., Patoka, R., et al. (2009). Induction of an Antiinflammatory Effect and Prevention of Cartilage Damage in Rat Knee Osteoarthritis by CF101 Treatment. Arthritis Rheum. 60 (10), 3061–3071. doi:10.1002/art.24817

PubMed Abstract | CrossRef Full Text | Google Scholar

Breedveld, F. C. (2004). Osteoarthritis-The Impact of a Serious Disease. Rheumatology (Oxford) 43 (Suppl 1), i4–i8. doi:10.1093/rheumatology/keh102

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., and Wang, X. (2020). miRDB: An Online Database for Prediction of Functional microRNA Targets. Nucleic Acids Res. 48 (D1), D127–D131. doi:10.1093/nar/gkz757

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, S., Fang, X.-Q., Zhang, J.-F., Ma, Y., Tang, X.-Z., Zhou, Z.-J., et al. (2016). Lycorine Protects Cartilage through Suppressing the Expression of Matrix Metalloprotenases in Rat Chondrocytes and in a Mouse Osteoarthritis Model. Mol. Med. Rep. 14 (4), 3389–3396. doi:10.3892/mmr.2016.5594

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Cui, S., Li, W., Jin, H., Fan, J., Sun, Y., et al. (2020). Ingenuity Pathway Analysis of Human Facet Joint Tissues: Insight into Facet Joint Osteoarthritis. Exp. Ther. Med. 19 (4), 2997–3008. doi:10.3892/etm.2020.8555

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, W.-S., Lin, H.-H., Ho, F.-M., Lai, C.-L., and Chao, C.-L. (2016). Risks of Acute Coronary Syndrome in Patients with Osteoarthritis: A Nationwide Population-Based Cohort Study. Clin. Rheumatol. 35 (11), 2807–2813. doi:10.1007/s10067-016-3391-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Cong, L., Zhu, Y., and Tu, G. (2017). A Bioinformatic Analysis of microRNAs Role in Osteoarthritis. Osteoarthr. Cartil. 25 (8), 1362–1371. doi:10.1016/j.joca.2017.03.012

CrossRef Full Text | Google Scholar

Cui, A., Li, H., Wang, D., Zhong, J., Chen, Y., and Lu, H. (2020). Global, Regional Prevalence, Incidence and Risk Factors of Knee Osteoarthritis in Population-Based Studies. EClinicalMedicine 29-30, 100587–30100587. doi:10.1016/j.eclinm.2020.100587

PubMed Abstract | CrossRef Full Text | Google Scholar

de Lange-Brokaar, B. J. E., Ioan-Facsinay, A., van Osch, G. J. V. M., Zuurmond, A.-M., Schoones, J., Toes, R. E. M., et al. (2012). Synovial Inflammation, Immune Cells and Their Cytokines in Osteoarthritis: a Review. Osteoarthr. Cartil. 20 (12), 1484–1499. doi:10.1016/j.joca.2012.08.027

CrossRef Full Text | Google Scholar

Eakin, G. S., Amodeo, K. L., and Kahlon, R. S. (2017). Arthritis and its Public Health Burden. Dela J. Public Health 3 (1), 36–44. doi:10.32481/djph.2017.03.006

CrossRef Full Text | Google Scholar

Fernandez-Moreno, M., Soto-Hermida, A., Pertega, S., Oreiro, N., Fernandez-Lopez, C., Rego-Perez, I., et al. (2011). Mitochondrial DNA (mtDNA) Haplogroups and Serum Levels of Anti-oxidant Enzymes in Patients with Osteoarthritis. BMC Musculoskelet. Disord. 12, 264. doi:10.1186/1471-2474-12-264

PubMed Abstract | CrossRef Full Text | Google Scholar

Guan, M., Zhu, Y., Liao, B., Tan, Q., Qi, H., Zhang, B., et al. (2020). Low‐intensity Pulsed Ultrasound Inhibits VEGFA Expression in Chondrocytes and Protects against Cartilage Degeneration in Experimental Osteoarthritis. FEBS Open Bio 10 (3), 434–443. doi:10.1002/2211-5463.12801

PubMed Abstract | CrossRef Full Text | Google Scholar

Hong, M., Tao, S., Zhang, L., Diao, L.-T., Huang, X., Huang, S., et al. (2020). RNA Sequencing: New Technologies and Applications in Cancer Research. J. Hematol. Oncol. 13 (1), 166. doi:10.1186/s13045-020-01005-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, K. (2020). Become Competent within One Day in Generating Boxplots and Violin Plots for a Novice without Prior R Experience. Methods Protoc. 3 (4), 64. doi:10.3390/mps3040064

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H.-Y., Lin, Y.-C. -D., Cui, S., Huang, Y., Tang, Y., Xu, J., et al. (2022). miRTarBase Update 2022: An Informative Resource for Experimentally Validated miRNA-Target Interactions. Nucleic Acids Res. 50 (D1), D222–D230. doi:10.1093/nar/gkab1079

PubMed Abstract | CrossRef Full Text | Google Scholar

Karagkouni, D., Paraskevopoulou, M. D., Tastsoglou, S., Skoufos, G., Karavangeli, A., Pierros, V., et al. (2020). DIANA-LncBase V3: Indexing Experimentally Supported miRNA Targets on Non-coding Transcripts. Nucleic Acids Res. 48 (D1), D101–D110. doi:10.1093/nar/gkz1036

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohl, M., Wiese, S., and Warscheid, B. (2011). Cytoscape: Software for Visualization and Analysis of Biological Networks. Methods Mol. Biol. 696, 291–303. doi:10.1007/978-1-60761-987-1_18

PubMed Abstract | CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R Package for Weighted Correlation Network Analysis. BMC Bioinform. 9 (1), 559–613. doi:10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Bai, B., Wang, J., Xu, Z., Yan, S., and Liu, G. (2017). Identification of Key mRNAs and microRNAs in the Pathogenesis and Progression of Osteoarthritis Using Microarray Analysis. Mol. Med. Rep. 16 (4), 5659–5666. doi:10.3892/mmr.2017.7251

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Zeng, Y., Si, H.-B., Tang, L., Xie, H.-Q., and Shen, B. (2022). Exosomes Derived from Human Urine-Derived Stem Cells Overexpressing miR-140-5p Alleviate Knee Osteoarthritis through Downregulation of VEGFA in a Rat Model. Am. J. Sports Med. 50 (4), 1088–1105. doi:10.1177/03635465221073991

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, E. B. P., Filiberti, A., Husain, S. A., and Humphrey, M. B. (2017). Immune Contributions to Osteoarthritis. Curr. Osteoporos. Rep. 15 (6), 593–600. doi:10.1007/s11914-017-0411-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, K., Shi, W., Xu, M., Liu, J., and Zhang, F. (2018). Genome-Wide Identification and Characterization of Long Non-coding RNA in Wheat Roots in Response to Ca²⁺ Channel Blocker. Front. Plant Sci. 9, 244. doi:10.3389/fpls.2018.00244

PubMed Abstract | CrossRef Full Text | Google Scholar

Metsios, G. S., Stavropoulos-Kalinoglou, A., Veldhuijzen van Zanten, J. J. C. S., Treharne, G. J., Panoulas, V. F., Douglas, K. M. J., et al. (2008). Rheumatoid Arthritis, Cardiovascular Disease and Physical Exercise: A Systematic Review. Rheumatol. Oxf. 47 (3), 239–248. doi:10.1093/rheumatology/kem260

CrossRef Full Text | Google Scholar

Muthukkumar, S., Nair, P., Sells, S. F., Maddiwar, N. G., Jacob, R. J., and Rangnekar, V. M. (1995). Role of EGR-1 in Thapsigargin-Inducible Apoptosis in the Melanoma Cell Line A375-C6. Mol. Cell Biol. 15 (11), 6262–6272. doi:10.1128/mcb.15.11.6262

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliviero, A., Della Porta, G., Peretti, G. M., and Maffulli, N. (2019). MicroRNA in Osteoarthritis: Physiopathology, Diagnosis and Therapeutic Challenge. Br. Med. Bull. 130 (1), 137–147. doi:10.1093/bmb/ldz015

PubMed Abstract | CrossRef Full Text | Google Scholar

Penkava, F., Velasco-Herrera, M. D. C., Young, M. D., Yager, N., Nwosu, L. N., Pratt, A. G., et al. (2020). Single-cell Sequencing Reveals Clonal Expansions of Pro-inflammatory Synovial CD8 T Cells Expressing Tissue-Homing Receptors in Psoriatic Arthritis. Nat. Commun. 11 (1), 4767. doi:10.1038/s41467-020-18513-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Poole, A. R., Ha, N., Bourdon, S., Sayre, E. C., Guermazi, A., and Cibere, J. (2016). Ability of a Urine Assay of Type II Collagen Cleavage by Collagenases to Detect Early Onset and Progression of Articular Cartilage Degeneration: Results from a Population-Based Cohort Study. J. Rheumatol. 43 (10), 1864–1870. doi:10.3899/jrheum.150917

PubMed Abstract | CrossRef Full Text | Google Scholar

Pundole, X., and Suarez-Almazor, M. E. (2020). Cancer and Rheumatoid Arthritis. Rheumatic Dis. Clin. N. Am. 46 (3), 445–462. doi:10.1016/j.rdc.2020.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Saberi, H. F., Siebuhr, A. S., Uitterlinden, A. G., Oei, E. H., Hofman, A., Karsdal, M. A., et al. (2016). Association between Biomarkers of Tissue Inflammation and Progression of Osteoarthritis: Evidence from the Rotterdam Study Cohort. Arthritis Res. Ther. 18, 81. doi:10.1186/s13075-016-0976-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Schieir, O., Tosevski, C., Glazier, R. H., Hogg-Johnson, S., and Badley, E. M. (2017). Incident Myocardial Infarction Associated with Major Types of Arthritis in the General Population: a Systematic Review and Meta-Analysis. Ann. Rheum. Dis. 76 (8), 1396–1404. doi:10.1136/annrheumdis-2016-210275

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, J., Abu-Amer, Y., O’Keefe, R. J., and McAlinden, A. (2017). Inflammation and Epigenetic Regulation in Osteoarthritis. Connect. Tissue Res. 58 (1), 49–63. doi:10.1080/03008207.2016.1208655

PubMed Abstract | CrossRef Full Text | Google Scholar

Skrzypa, M., Szala, D., Gablo, N., Czech, J., Pajak, J., Kopanska, M., et al. (2019). miRNA-146a-5p is Upregulated in Serum and Cartilage Samples of Patients with Osteoarthritis. Pol. Przegl Chir. 91 (3), 1–5. doi:10.5604/01.3001.0013.0135

PubMed Abstract | CrossRef Full Text | Google Scholar

Smyth, G. K., Michaud, J., and Scott, H. S. (2005). Use of Within-Array Replicate Spots For Assessing Differential Expression In Microarray Experiments Bioinformatics 21 (9), 2067–2075. doi:10.1093/bioinformatics/bti270

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. U.S.A. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Tavallaee, G., Rockel, J. S., Lively, S., and Kapoor, M. (2020). MicroRNAs in Synovial Pathology Associated with Osteoarthritis. Front. Med. (Lausanne) 7, 376. doi:10.3389/fmed.2020.00376

PubMed Abstract | CrossRef Full Text | Google Scholar

Tian, Z., He, W., Tang, J., Liao, X., Yang, Q., Wu, Y., et al. (2020). Identification of Important Modules and Biomarkers in Breast Cancer Based on WGCNA. Onco Targets Ther. 13, 6805–6817. doi:10.2147/OTT.S258439

PubMed Abstract | CrossRef Full Text | Google Scholar

Ungethuem, U., Haeupl, T., Witt, H., Koczan, D., Krenn, V., Huber, H., et al. (2010). Molecular Signatures and New Candidates to Target the Pathogenesis of Rheumatoid Arthritis. Physiol. Genomics 42A (4), 267–282. doi:10.1152/physiolgenomics.00004.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Raemdonck, K., Umar, S., Palasiewicz, K., Volkov, S., Volin, M. V., Arami, S., et al. (2020). CCL21/CCR7 Signaling in Macrophages Promotes Joint Inflammation and Th17-Mediated Osteoclast Formation in Rheumatoid Arthritis. Cell. Mol. Life Sci. 77 (7), 1387–1399. doi:10.1007/s00018-019-03235-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, C., Zhang, F., Yao, H., Li, H., and Tuan, R. S. (2021). Histone Modifications and Chondrocyte Fate: Regulation and Therapeutic Implications. Front. Cell Dev. Biol. 9, 626708. doi:10.3389/fcell.2021.626708

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Yu, Y., Huang, Y., Zhu, M., Chen, R., Liao, Z., et al. (2020). Identification of Potential Diagnostic Gene Biomarkers in Patients with Osteoarthritis. Sci. Rep. 10 (1), 13591. doi:10.1038/s41598-020-70596-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Watt, F. E. (2018). Osteoarthritis Biomarkers: Year in Review. Osteoarthr. Cartil. 26 (3), 312–318. doi:10.1016/j.joca.2017.10.016

CrossRef Full Text | Google Scholar

Woetzel, D., Huber, R., Kupfer, P., Pohlers, D., Pfaff, M., Driesch, D., et al. (2014). Identification of Rheumatoid Arthritis and Osteoarthritis Patients by Transcriptome-Based Rule Set Generation. Arthritis Res. Ther. 16 (2), R84. doi:10.1186/ar4526

PubMed Abstract | CrossRef Full Text | Google Scholar

World Medical Association (2013). World Medical Association Declaration of Helsinki: Ethical Principles for Medical Research Involving Human Subjects. JAMA 310 (20), 2191–2194. doi:10.1001/jama.2013.281053

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, H., Wang, C., Song, H., Xu, Y., and Ji, G. (2019). RNA-seq Profiling of Circular RNAs in Human Colorectal Cancer Liver Metastasis and the Potential Biomarkers. Mol. Cancer 18 (1), 8. doi:10.1186/s12943-018-0932-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: An R Package for Comparing Biological Themes Among Gene Clusters. Omics J. Integr. Biol. 16 (5), 284–287. doi:10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, W.-H., Xie, Q.-Q., Wang, K.-P., Shen, W., Feng, X.-F., Liu, Z., et al. (2021). Screening of Osteoarthritis Diagnostic Markers Based on Immune-Related Genes and Immune Infiltration. Sci. Rep. 11 (1), 7032. doi:10.1038/s41598-021-86319-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, F., He, Q., Tsang, W. P., Garvey, W. T., Chan, W. Y., and Wan, C. (2014). Insulin Exerts Direct, IGF-1 Independent Actions in Growth Plate Chondrocytes. Bone Res. 2, 14012. doi:10.1038/boneres.2014.12

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, R., Guo, H., Yang, X., Li, Z., Zhang, D., Li, B., et al. (2019). Potential Candidate Biomarkers Associated with Osteoarthritis: Evidence from a Comprehensive Network and Pathway Analysis. J. Cell. Physiology 234 (10), 17433–17443. doi:10.1002/jcp.28365

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, H., Liu, R., Sun, L., Guo, W., Ji, X., and Hu, X. (2021). Comprehensive Analysis of Gene Expression Changes and Validation in Hepatocellular Carcinoma. Onco Targets Ther. 14, 1021–1031. doi:10.2147/ott.s294500

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., and Xu, J. (2018). Synovial Fluid-Derived Exosomal lncRNA PCGEM1 as Biomarker for the Different Stages of Osteoarthritis. Int. Orthop. (SICOT) 42 (12), 2865–2872. doi:10.1007/s00264-018-4093-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Lv, J., Zhang, H., Xie, J., Dai, H., and Zhang, X. (2021). Gene Expression Profiles Analyzed Using Integrating RNA Sequencing, and Microarray Reveals Increased Inflammatory Response, Proliferation, and Osteoclastogenesis in Pigmented Villonodular Synovitis. Front. Immunol. 12, 665442. doi:10.3389/fimmu.2021.665442

PubMed Abstract | CrossRef Full Text | Google Scholar

Zupan, J., Vrtačnik, P., Cör, A., Haring, G., Weryha, G., Visvikis-Siest, S., et al. (2018). VEGF-A is Associated with Early Degenerative Changes in Cartilage and Subchondral Bone. Growth factors. 36 (5-6), 263–273. doi:10.1080/08977194.2019.1570926

PubMed Abstract | CrossRef Full Text | Google Scholar