{'fraction_non_alphanumeric': 0.036756794279962514, 'fraction_numerical': 0.024435111589323523, 'mean_word_length': 5.2730241672109734, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 8, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Single-Cell Profiles of Age-Related Osteoarthritis Uncover Underlying Heterogeneity Associated With Disease Progression

Objective: Osteoarthritis (OA) is the most common chronic degenerative joint disease, which represents the leading cause of age-related disability. Here, this study aimed to depict the intercellular heterogeneity of OA synovial tissues. Methods: Single-cell RNA sequencing (scRNA-seq) data were preprocessed and quality controlled by the Seurat package. Cell cluster was presented and cell types were annotated based on the mRNA expression of corresponding marker genes by the SingleR package. Cell-cell communication was assessed among different cell types. After integrating the GSE55235 and GSE55457 datasets, differentially expressed genes were identified between OA and normal synovial tissues. Then, differentially expressed marker genes were overlapped and their biological functions were analyzed. Results: Totally, five immune cell subpopulations were annotated in OA synovial tissues including macrophages, dendritic cells, T cells, monocytes and B cells. Pseudo-time analysis revealed the underlying evolution process in the inflammatory microenvironment of OA synovial tissue. There was close crosstalk between five cell types according to the ligand-receptor network. The genetic heterogeneity was investigated between OA and normal synovial tissues. Furthermore, functional annotation analysis showed the intercellular heterogeneity across immune cells in OA synovial tissues. Conclusion: This study offered insights into the heterogeneity of OA, which provided in-depth understanding of the transcriptomic diversities within synovial tissue. This transcriptional heterogeneity may improve our understanding on OA pathogenesis and provide potential molecular therapeutic targets for OA.


INTRODUCTION
Osteoarthritis (OA) represents the most frequent ageassociated chronic degenerative joint disease (Maumus et al., 2020). Aging is the key driving force in OA (Deng et al., 2019). The incidence of OA in the middle-aged population is 40-80%, and the disability rate is >50% (Zhang et al., 2020b). The diverse aetiology of OA is usually the result of many overlapping factors. The pathogenesis of OA involves the entire joint, composed of articular cartilage, synovial membrane as well as subchondral bone (Griffin and Scanzello, 2019). Traditional therapeutic strategies such as non-steroidal anti-inflammatory drugs only ease pain symptoms. Joint replacement may effectively treat OA patients at the end stage, with relatively high surgery risk and economic cost (Zhou et al., 2020). At current, there is still a lack of FDA-approved disease-modifying OA drugs that may remiss the pain and restrain joint degradation (Grandi and Bhutani, 2020).
Synovitis is common in OA and may augment cartilage damage, which has been a therapeutic target of OA (Labinsky et al., 2020). Synovitis is characterized by synovial tissue hyperplasia, inner macrophage aggregation as well as cell secretion disfunction ). The immune system that is involved in the initiation and progression of OA is a key element in the pathogenesis of this disease . The inflammatory microenvironment of OA synovial tissues consists of immune cells and inflammatory mediators such as inflammatory cytokines TNFα, and anti-inflammatory cytokines (IL-4, IL-10, and IL-13) (Wojdasiewicz et al., 2014). The pathophysiological processes that occur in the OA joint are mostly mediated by inflammatory mediators secreted from immune cells. Determining the characteristics of immune cells and mediators in OA synovium may reveal critical features of OA pathogenesis (Zheng et al., 2021). So far, it remains indistinct on the roles of inflammatory mediators in the context of OA synovium.
Single-cell RNA sequencing (scRNA-seq) has been applied for uncovering molecular programs and lineage progression of OA (Zhang et al., 2020c). For instance, Zhang et al. (2020c) identified the biomarkers and differentiation of chondrocytes in human OA using scRNA-seq analyses. Moreover, Ji et al. (2019) adopted scRNA-seq analyses to reveal chondrocyte taxonomy and mechanisms of human OA cartilage regeneration. Despite this, in-depth exploration should be conducted for clarifying the pathogenesis of OA. In this study, we identified differentially expressed genes (DEGs) in OA synovium and found that these DEGs mainly participated in inflammatory pathways based on transcriptome data. Moreover, we revealed immune cell subpopulations and their intercellular heterogeneity in OA synovial tissues by single-cell RNA sequencing (scRNA-seq) profiles. These findings might deepen our understanding on the transcriptomic diversities within synovial tissue.

Differential Expression Analyses
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.

Functional Annotation Analysis
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.

CIBERSORT Analysis
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.

ScRNA-Seq and Preprocessing
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.

Cell Cluster
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).

Pseudo-Time Analysis
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.

Ligand-Receptor Network Analysis
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. Ligandreceptor relationships between cells were then identified. Finally, ligand-receptor network was conducted via Cytoscape software (version 3.7.1) (Shannon et al., 2003).

Identification of Differentially Expressed Marker Genes
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.

Cell Culture and Establishment of OA Cellular Models
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 2 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.

Flow Cytometry
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).

5-Ethynyl-29-Deoxyuridine (EdU) Staining
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).

Statistical Analyses
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).

Screening DEGs in OA Synovial Tissues
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 Figure 1E; Table 1).

Functional Annotation Analysis of Common DEGs
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).

Landscape in Immune Infiltration Between OA and Healthy Synovial Tissues
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).

Differences in Immune Infiltration Between OA and Healthy Synovial Tissues
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).

Quality Control of scRNA-Seq of Osteoarthritic Synoviocytes
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).

scRNA-Seq Uncovers Six Cell Types for OA Synoviocytes
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 Frontiers in Molecular Biosciences | www.frontiersin.org January 2022 | Volume 8 | Article 748360 8 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.

Identification of Five Immune Cell Subpopulations and Associated Marker Genes for OA Synoviocytes
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,  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.

Pseudo-Time Analysis of Osteoarthritic Synoviocytes
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.

Ligand-Receptor Network
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.

Biological Functions of Differentially Expressed Marker Genes
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 upregulated genes and 91 down-regulated genes in OA compared with normal synovial tissues ( Figure 10C 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.

Silencing NR4A1 and NR4A2 Relieves IL-1β-Induced Chondrocyte Damage
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 upregulations 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.

DISCUSSION
OA represents an age-related arthritis and the major chronic disability-related cause among elderly populations (Zhang et al., 2020a). OA is characterized by synovial inflammation, confirmed by magnetic resonance imaging and ultrasound (Conaghan et al., 2019). Here, based on microarray and scRNA-seq profiles, we comprehensively expounded immune cells and inflammatory mediators in OA synovium, which characterized the key features of OA pathogenesis. More importantly, our study depicted the intercellular heterogeneity of OA synovial tissues.
With the advance in RNA-seq and microarray techniques, more and more studies have revealed the molecular features of OA. Previous findings on gene expression profiles of OA synovial tissues have proposed several key genes as well as their enriched signaling pathways (Wang et al., 2020). Herein, this study identified 57 DEGs in OA compared to healthy synovial tissues. These genes were involved in several inflammatory pathways such as TNF, IL-17, NF-κB, Toll-like receptor and Th17 cell differentiation, indicating that they might be related to OA progression. Especially, as reported in a previous study, OArelated DEGs were also distinctly enriched in TNF pathway . IL-17 regulates response to damage in OA by immunologically inducing senescence (Faust et al., 2020). NF-κB pathway affects synovial inflammation, which leads to the damage of joint cartilage (Choi et al., 2019). Moreover, Tolllike receptor may contribute to pathological alterations in OA joint tissue specimens (Miller et al., 2019). IL-17 mediates apoptosis in rheumatoid arthritis synoviocytes by activating autophagy (Kim et al., 2017). Therefore, alterations in gene expression may participate in inflammatory pathways of OA.
The inflammatory microenvironment of the OA synovium consists of immune cells and various inflammatory cytokines. The CIBERSORT algorithm confirmed that there were relatively high infiltration levels of macrophages in OA synovial tissues. Consistently, lower infiltration levels of dendritic cells activated, T cells CD4 memory resting, neutrophils and mast cells activated were detected in OA than healthy synovial tissues (Cai et al., 2020;Chen et al., 2020). On the contrary, there were higher infiltration levels of macrophages M0, mast cells resting and monocytes in OA compared to control synovial tissue specimens. These data were indicative that dysregulated immune cells in synovium might play key roles in OA pathogenesis. To further understand the roles of immune cells in the progression of OA, our scRNA-seq data demonstrated that OA may involve five immune cells including macrophages, dendritic cells, T cells, monocytes and B cells as well as associated effector cytokines and cytokine receptors in synoviocytes. The ligand-receptor network revealed the close crosstalk between cells. Among them, macrophages were most abundant immune cells. Consistently, Hsueh et al. (2021) found that macrophages were predominant in the OA synovial tissues, occupying 53% of total synoviocytes. Macrophages could product inflammatory cytokines including TNF, IL1B, IL18, TGFB1, and IL8. It has been confirmed that macrophages and macrophage-secreted mediators force the inflammatory and destructive responses in the OA synovial tissues (Bondeson et al., 2006). Therefore, macrophages may participate in the pathogenesis of OA. Inflammatory cytokines in the synovial membrane are also involved in the initiation and progression of OA . For example, TNF-α induces the production of matrix metalloproteinases like MMP-13, ADAMTS-5 and ADAMTS-7 and activates the NF-κB pathway that can facilitate the proinflammatory functions of TNF-α (Zhao et al., 2019). Above data demonstrated that the immune cells, cytokines and receptors in synovial tissues might play critical roles in OA progression. Our in vitro experiments confirmed that silencing differentially expressed marker genes NR4A1 and NR4A2 relieved IL-1β-induced chondrocyte damage in vitro experiments. However, a few limitations should be pointed out. First, more scRNA-seq data should be verified for the intercellular heterogeneity of OA synovial tissues. Second, in vivo experiments should be carried out to verify the biological functions of NR4A1 and NR4A2 in OA progression.

CONCLUSION
This study uncovered the transcriptomic diversities within OA synovial tissue. Our findings could deepen the understanding on the pathophysiology of OA and facilitate the progression of novel drugs against therapeutic targets.

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

AUTHOR CONTRIBUTIONS
WS and MM conceived and designed the study. WL, YC, and GZ conducted most of the experiments and data analysis, and wrote the manuscript. SY and TY participated in collecting data and helped to draft the manuscript. All authors reviewed and approved the manuscript.