Gene expression profiles of GSE55235 and GSE55457 were retrieved from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/gds/) using the GEOquery package (version 3.6.1; https://www.r-project.org/). The GSE55235 dataset contained microarray expression profiling of 10 synovial tissues from healthy joints and 10 synovial tissues from osteoarthritic joints (Woetzel et al., 2014). The GSE55457 dataset contained expression profiling of 10 synovial tissue specimens from normal joints and 10 synovial tissue specimens from osteoarthritic joints (Woetzel et al., 2014). Above datasets were based on GPL96 platform, [HG-U133A] Affymetrix Human Genome U133A Array.

Raw CEL files were read by affy package (version 1.0) (Gautier et al., 2004). Microarray expression profiling was normalized with Robust MultiChip Analysis (RMA) method (Irizarry et al., 2003) and standardized by quantile. The expression value was then log2 converted. The Linear Models for Microarray Data (limma; version 1.0) package was applied for screening DEGs between OA and healthy synovial tissues with student’s t-test in the GSE55235 and GSE55457 datasets (Ritchie et al., 2015). False discovery rate (FDR) was determined for multiple comparisons with Benjamini and Hochberg method. Genes with FDR<0.05 and |fold-change|>1.5 were defined as DEGs. Hierarchical clustering analyses of DEGs were presented with the gplots package. DEGs in the two datasets were intersected to obtain common DEGs for OA.

The biological functions and pathways of common DEGs were interpreted via the clusterProfiler package (version 1.0; http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) (Yu et al., 2012). Biological processes of gene ontology (GO) were functionally annotated for these DEGs. For understanding signaling pathways enriched by these DEGs, Kyoto encyclopedia of genes and genomes (KEGG) pathway enrichment analysis was employed. Terms with FDR<0.05 were considered significant enrichment. The top ten biological processes and the top 30 KEGG pathways were visualized.

The gene expression profiles of the two datasets were integrated into one dataset. The CIBERSORT algorithm (version 1.0; http://cibersort.stanford.edu/) was applied to characterize immune cell compositions in each specimen based on normalized gene expression profiling (Newman et al., 2015). The specimens with p < 0.05 were screened and the proportions of 22 kinds of immune cells were determined, including B cells naïve, B cells memory, plasma cells, T cells CD8, T cells CD4 naïve, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells regulatory (Tregs), T cells gamma delta, NK cells resting, NK cells activated, monocytes, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, eosinophils and neutrophils. Principal component analysis (PCA) was presented for determining the differences in the proportions of immune cells between OA and healthy synovial tissue specimens. The infiltration levels of immune cells between groups were compared using the vioplot package.

Single-cell transcriptomic data of 10,640 synoviocytes from 3 osteoarthritic synovial membrane samples were obtained from the GSE152805 dataset on the GPL20301 Illumina HiSeq 4000 (Homo sapiens) (Chou et al., 2020). The Seurat package (version 3.0; http://satijalab.org/seurat/) was employed for quality control (Butler et al., 2018). Cells with 200–5000 feature genes and the percent of transcripts of mitochondrial genes <10% were retained for further analysis. Expression data were normalized with the LogNormalize method. Briefly, based on the total expression levels, the gene expression value in each cell was normalized and then multiplied by a scale factor = 10,000, followed by log2 conversion. Gene variances were calculated via the FindVariableFeatures function. Each dataset returned the top 2000 highly variable genes, which were used for downstream analysis. Meanwhile, the top 10 genes with the highest variances were selected and visualized. Then, the data were scaled utilizing the ScaleData function, as follows: the expression of each gene was shifted so that the average expression between cells was 0 and the expression of each gene was scaled so that the difference between cells was 1. After scaling the data, PCA was presented for the 2000 highly variable genes with the VizDimReduction and DimPlot functions. The number of principal components (PCs) was determined with the elbow plot method.

In PCA, a KNN graph was constructed based on Euclidean distance. The FindNeighbors function was applied to find their edge weights between any 2 cells. Cell cluster was analyzed based on the FindClusters function. Non-linear dimensionality reduction technology t-distributed statistical neighbor embedding (tSNE) was employed to cluster similar cells together in a low-dimensional space. Markers were identified for each cluster using the FindAllMarkers function, which were visualized with the DoHetmap, VlinPlot and FeaturePlot functions. Based on the markers of each cluster, cell type was annotated via the Cell Type Annotation of Single Cells (SingleR; version 1.0) (Aran et al., 2019).

Monocle algorithm (version 2.0) was applied for pseudo-time analysis (Qiu et al., 2017). Genes that were expressed in >5% of the cells were selected. The t-SNE was presented with reduce Dimension function and cells were clustered with cluster Cells function. Genes that were differentially expressed (p-value < 0.05) between clusters were selected as candidate genes with differential Gene Test function. Through the reduce Dimension function, based on the above candidate genes, dimensionality reduction analysis of the cells was caried out with DDRTree method. Then, cells were sorted and visualized according to single-cell trajectories using order Cells function.

By CellPhoneDB (version 2.0) project (Efremova et al., 2020), the cell-cell communication was estimated based on the expression of receptors and ligands corresponding to each cell type. Ligand-receptor relationships between cells were then identified. Finally, ligand-receptor network was conducted via Cytoscape software (version 3.7.1) (Shannon et al., 2003).

The GSE55235 and GSE55457 datasets were integrated and batch effects were corrected. Differential expression analysis was carried out between OA and normal synovial tissues. Genes with adjusted p < 0.05 and |fold-change|>1.5 were considered differentially expressed. By overlapping marker genes in each cell type, differentially expressed marker genes were identified. Functional annotation analysis of above marker genes was performed.

Human chondrocytes C28/I2 (ATCC, United States) were grown in DMEM (Sigma, United States) as well as were incubated in a humidified environment with 10% fetal bovine serum containing 5% CO₂ at 37°C. Chondrocytes were planted into a six-well plate and exposed to 10 ng/ml IL-1β (Sigma, United States) lasting 24 h to establish OA cellular models.

Chondrocytes were lysed with RIPA buffer (pH 8.0) and extracted protein was quantified with BCA Protein Quantitation Kit. Thereafter, protein was isolated on SDS-PAGE gels as well as transferred onto PVDF membrane. The membrane was sealed lasting 1 h utilizing PBST plus 5% BSA and incubated by specific primary antibodies at 4°C overnight, containing NR4A1 (1:500; #25851-1-AP; Proteintech, Wuhan, China), NR4A2 (1:500; #66878-1-Ig; Proteintech), GAPDH (1:1000; #60004-1-Ig; Proteintech), MMP3 (1:500; #17873-1-AP; Proteintech), MMP13 (1:1000; #18165-1-AP; Proteintech) and Collagen II (1:800; #28459-1-AP; Proteintech). Afterwards, the membrane was incubated by HRP-conjugated goat anti-rabbit secondary antibody (1:2000; #SA00001-2; Proteintech) lasting 1 h. Protein bands were visualized utilizing ECL kits.

Small interfering RNAs (siRNAs) of NR4A1 (si-NR4A1) and NR4A2 (si-NR4A2) were synthesized by GenePharma (Shanghai, China). IL-1β-induced chondrocytes were seeded onto a six-well plate (1×10⁶/well). When cell confluence reached 80–90%. In line with the manufacturer’s protocol, cell transfection was presented utilizing Lipofectamine 2000 Reagent (Invitrogen, United States). After 48 h, NR4A1 and NR4A2 expressions were verified with western blotting.

Chondrocyte apoptosis was determined utilizing Annexin V-FITC Apoptosis Detection kits (BD, United States) in accordance with the manufacturer’s protocol. Chondrocytes were collected and washed twice by PBS. Thereafter, chondrocytes were suspended through adding Annexin V, followed by incubation at 4°C lasting 15 min in the dark. Chondrocytes were incubated by propidium iodide lasting 5 min in the dark. Stained chondrocytes were tested utilizing flow cytometer (FACSCalibur; BD, United States).

EdU cellular proliferation experiment was conducted. In brief, pre-warmed 2× EdU working solution (RiboBio, China) was added to the cell plate to make EdU 10 Μm concentration, followed by incubation lasting 2 h. Thereafter, chondrocytes were fixed by 4% paraformaldehyde as well as permeated by 0.3% Triton X-100 solution at room temperature lasting 15 min. DAPI (Sigma, United States) was then added to each well as well as incubated lasting 10 min. Images were captured under a Nikon A1Si Laser Scanning Confocal Microscope (Nikon, Japan).

Data were displayed as the mean ± SD, and statistical analyses were carried out using Student’s t test or one-way ANOVA for multiple comparisons. p < 0.05 was statistically significant. All statistical analyses were implemented utilizing R software (version 3.6.1; https://www.r-project.org) and GraphPad Prism (version 8.0.1).

The gene expression profiles of OA and healthy synovial tissue specimens were obtained from two microarray datasets GSE55235 and GSE55457. With the criteria of FDR<0.05 and |fold-change|>1.5, 148 genes were up-regulated and 127 genes were down-regulated in OA compared to healthy synovial tissues in the GSE55235 dataset (Figures 1A,B; Supplementary Table S1). Meanwhile, there were 98 up- and 91 down-regulated genes in OA synovial tissues in the GSE55457 dataset (Figures 1C,D; Supplementary Table S2). To obtain genes closely related to OA, 57 common DEGs were taken intersection in the two datasets, as follows: EMP1, PTN, ITGA5, IL10RA, FOSL2, EPHA3, ZFP36, STC1, TRIL, CHD1, HTR2B, ACOT13, SFPQ, CCNL1, NAP1L3, B3GALNT1, GLT8D2, CASP1, GTF2H5, EML1, PTGS2, YBX3, HPGDS, THBS4, SMCO4, NR4A1, FKBP5, TLR7, DDX3Y, NFKBIA, JUNB, PTGS1, MYL6B, SLC2A3, SAMSN1, XIST, LONRF1, TBL1XR1, VEGFA, SLC5A3, EFNB2, MYC, JUN, PNMAL1, PSMB10, KCNN4, PRSS23, NTAN1, TM6SF1, MAFF, TNFAIP3, MCL1, TMEM230, HLA-DQB1, GADD45B, IL11RA and RLF (Figure 1E; Table 1).

To uncover biological functions and pathways of these 57 common DEGs, functional enrichment analyses were carried out. GO enrichment results showed that these DEGs were involved in regulating epithelial cell proliferation and response to mechanical stimulus, lipopolysaccharide and molecule of bacterial origin (Figure 2A; Table 2). Moreover, several key pathways were distinctly enriched by these DEGs such as TNF, IL-17, NF-kappa B, apoptosis, osteoclast differentiation, PI3K-Akt, JAK-STAT, Toll-like receptor and Th17 cell differentiation pathways (Figure 2B; Table 2).

We employed the CIBERSORT algorithm to evaluate the immune cell infiltrations in OA and healthy synovial tissues from the integrated GSE55235 and GSE55457 datasets. Our PCA results demonstrated that the infiltration levels of immune cell compositions from OA and healthy synovial tissues exhibited significant group-bias cluster as well as personal disparity (Figure 3A). As shown in Figure 3B, we summarized the relative percent of different immune cells in each sample. Heatmap depicted that the infiltration levels of macrophages M2 were relatively high both in OA and healthy synovial tissue samples (Figure 3C). We further found that there were distinct correlations between immune cells (Figure 3D). Especially, T cells CD4 naïve were highly correlated to neutrophils (r = 0.91).

This study compared the differences in infiltration levels of immune cells between OA and healthy synovial tissues from the GSE55235 and GSE55457 datasets with the CIBERSORT method (Figure 4A). Lower infiltration levels of dendritic cells activated (p = 0.035), T cells CD4 memory resting (p = 0.034), neutrophils (p = 0.001) and mast cells activated (p = 1.156e-04) were found in OA synovial tissues compared to healthy synovial tissues (Figures 4B–E). Oppositely, there were higher levels of macrophages M0 (p = 0.015), mast cells resting (p = 2.105e-06) and monocytes (p = 0.049) in OA than control synovial tissue specimens (Figures 4F–H).

This study analyzed the scRNA-seq data of osteoarthritic synoviocytes from the GSE152805 dataset. We detected the genes in each cell and removed low-quality cells or empty droplets with very few genes <500 and cells with doublets or multiplets with high gene counts >2000. In addition, we filtered out the cells with percent of the mitochondrial genome >5%. After filtration, we visualized the number of feature genes and the count of genes as well as the percent of mitochondrial genes in OA synoviocytes (Figure 5A). Also, there was very low correlation between the count of genes and the percent of mitochondrial genes (Figure 5B). The count of genes was highly positively correlated to the number of feature genes. After filtration, the data were normalized with the LogNormalize method and the top 2000 highly variable genes were identified in OA synoviocytes for downstream analysis (Figure 5C). The top ten genes were as follows: TPSAB1, PENK, AMTN, IL1B, CCL4, CCL20, SELE, CXCL14, CCL3, and CSN1S1. After dimensionality reduction analysis, we visualized the highly variable genes and OA synoviocytes in the top two PCs (Figures 5D,E). With the elbow plot method, the top 30 PCs were determined for further analysis (Figure 5F).

Totally, 10,154 OA synoviocytes were clustered into 14 cell populations via the t-SNE method (Figure 6A). As shown in Figure 6B, we visualized the top 10 differentially expressed marker genes for each cell cluster. With the SingleR package, five cell types were annotated, as follows (from most to least): chondrocytes, macrophages, monocytes, fibroblasts, adipocytes and T cells (Figure 6C). Also, the cell cycle distribution of each cell cluster was shown (Figure 6D). Secreted inflammatory cytokines and their receptors are considered as key mediators for OA progression (Kapoor et al., 2011). Here, we analyzed the expression distribution of cytokine receptors in cell types. As a result, TGFBR1 was mainly expressed in macrophages and TNFRSF11B was mainly expressed in chondrocytes (Figures 6E–G). Meanwhile, TNFRSF1B was primarily distributed in monocytes, followed by T cells and macrophages. IL1R1 was mainly expressed in fibroblasts and chondrocytes, while IL10RA was primarily distributed in T cells, followed by monocytes and macrophages. Also, TNFRSF1A was expressed in fibroblasts, adipocytes, chondrocytes, macrophages and monocytes. This study further visualized the distribution of secreted inflammatory cytokines in each cell type (Figures 6H–J). TNF and IL-18 were expressed in macrophages and monocytes as well as IL1A and IL10 were only expressed in monocytes. IL-6 was primarily expressed in fibroblasts and IL1B was mainly distributed in monocytes, followed by macrophages. TGFB1 was mainly expressed in T cells, followed by monocytes, macrophages, adipocytes, fibroblasts and chondrocytes. Meanwhile, IL8 was primarily distributed in macrophages, followed by monocytes, fibroblasts, T cells, adipocytes and chondrocytes.

Next, we further annotated the immune cell subpopulations of macrophages, monocytes and T cells across 1353 OA synoviocytes via the SingleR package. As a result, five cell subpopulations were obtained, including macrophages, dendritic cells, T cells, monocytes and B cells (from most to least; Figure 7A). The top 20 marker genes for each cell subpopulation were visualized, as shown in Figure 7B. The expression distributions of cytokine receptors were detected in each cell subpopulation, as follows: TGFBR1 (monocytes), TNFRSF11B (monocytes and B cells), TNFRSF1B (dendritic cells, T cells, monocytes and macrophages), IL1R1 (monocytes), IL10RA (T cells, monocytes, dendritic cells, macrophages, B cells and monocytes) and TNFRSF1A (monocytes, macrophages and dendritic cells; Figures 7C,D). The expression of inflammatory cytokines was shown in each cell subpopulation, as follows: TNF (dendritic cells, monocytes and macrophages), IL1A (dendritic cells and monocytes), IL6 (monocytes and dendritic cells), IL1B (dendritic cells, monocytes, macrophages and B cells), IL18 (macrophages, dendritic cells and monocytes), TGFB1 (T cells, dendritic cells, monocytes, macrophages and B cells), IL10 (dendritic cells and monocytes) and IL8 (macrophages, dendritic cells, monocytes, B cells and T cells; Figures 7E,F). Also, we found that CD79A was mainly expressed in B cells, CD3E was mainly expressed in T cells and ITGAX was mainly expressed in dendritic cells (Figure 7G). Meanwhile, ITGAM was mainly expressed in macrophages, dendritic cells and monocytes.

Monocle 2 algorithm was used to reshape the change process of cells over time by constructing the trajectory of changes between cells. Figure 8A showed the distribution of osteoarthritic synoviocytes in 7 different states. As shown in Figure 8B, we observed the heterogeneity in the mRNA expression of inflammatory cytokines (IFNG, IL10, IL18, IL1A, IL1B, IL6, IL8, and TNF) and cytokine receptor (TNFRSF1A) among states. The distribution of osteoarthritic synoviocytes according to pseudo-time value was shown in Figure 8C. The scatter plot demonstrated the dynamic expression of above genes over the pseudo-time value (Figure 8D). As the pseudo-time value increased, the mRNA expression of IFNG gradually decreased while the expression of IL10, IL18, IL1A, IL1B, IL6, IL8, TNF, and TNFRSF1A gradually increased. Above findings revealed the developmental characteristics of osteoarthritic synoviocytes.

By CellPhoneDB v2.0, this study identified receptor-ligand relationships between five cell types. Figure 9 showed that the close interactions between cell types according to receptor-ligand interactions. The arrow pointed to the recipient cell. The width of the line represented the number of ligand-receptor relationships.

We integrated the GSE55235 and GSE55457 datasets and removed batch effects (Figures 10A,B). With the criteria of FDR<0.05 and |fold-change|>1.5, we identified 121 up-regulated genes and 91 down-regulated genes in OA compared with normal synovial tissues (Figure 10C). After overlapping marker genes of each cell type, we finally identified 6 differentially expressed marker genes in B cells (ASNS, NDUFA3, SPARC, GCHFR, BHLHE41, and BGN), 22 differentially expressed marker genes in dendritic cells (CXCL3, FOSL2, MAP3K8, FCER1A, PTGS2, OLR1, KDM6B, PFKFB3, NFKBIA, GCH1, SAMSN1, HLA-DPB1, VEGFA, NR4A3, KCNN4, PMAIP1, MAFF, TNFAIP3, FOSL1, HLA-DQB1, CD58, and GLRX), 27 differentially expressed marker genes in macrophages (VAMP8, NR4A2, DNAJB1, HNMT, TREM2, NR4A1, NTAN1, GRN, TMEM230, MNDA, KLF4, EMP1, CRIP1, SLC7A7, CEBPD, GDE1, NCF1, HPGDS, CD83, FKBP5, KLF10, FOLR2, SGMS1, JUN, TMEM176B, GADD45B, and FRMD4B), 35 differentially expressed marker genes in monocytes (FERMT2, DPT, FAP, PCOLCE, PTGDS, GLT8D2, ANKH, TMCO3, THBS4, DDX3Y, JUNB, TNFRSF11B, PTPRG, PRSS23, IFI27, VCAN, GSN, TGFBR3, ANGPTL2, SSPN, THY1, PPIC, SPARC, TPPP3, FZD1, MYL6B, MYC, MT1X, CRISPLD2, BGN, COL1A2, IL1R1, LTBP3, APOC1, and GADD45B), and 9 differentially expressed marker genes in T cells (CD37, SLC2A3, SAMSN1, CD52, ZFP36, ELF1, PRRC2C, MT1X, and SCAF11). Functional annotation analysis was used for investigating their biological functions. We observed that differentially expressed marker genes in B cells were mainly involved in modulating the molecular functions of extracellular matrix (Figure 10D). Differentially expressed marker genes in dendritic cells markedly participated in immune-related processes and pathways (Figures 10E,F). In Figure 10G, differentially expressed marker genes in macrophages were mainly involved in modulating inflammatory response and leukocyte differentiation. Differentially expressed marker genes in monocytes primarily participated in osteoclast differentiation (Figures 10H,I). In Figure 10J, differentially expressed marker genes in T cells were mainly involved in regulating cytoplasmic stress granule.

We further verified the biological functions of differentially expressed marker genes NR4A1 and NR4A2 in OA progression. IL-1β-induced OA chondrocyte models were constructed and our results confirmed the remarkable up-regulations of NR4A1 and NR4A2 in IL-1β-induced chondrocytes (Figures 11A–C). The expressions of NR4A1 and NR4A2 were then silenced by specific siRNAs in IL-1β-induced chondrocytes (Figures 11D–F). Flow cytometric analyses demonstrated that knockdown of NR4A1 and NR4A2 relieved IL-1β-induced chondrocyte apoptosis (Figures 11G–I). Additionally, silencing NR4A1 and NR4A2 enhanced proliferative capacity of IL-1β-induced chondrocyte apoptosis (Figures 11J–L). We also investigated that knockdown of NR4A1 and NR4A2 reduced the expressions of MMP3 and MMP3 as well as increased the expression of Collagen II in IL-1β-induced chondrocyte apoptosis (Figures 11M–T). Overall, silencing NR4A1 and NR4A2 relieved IL-1β-induced chondrocyte damage.