Single-Cell Profiles of Age-Related Osteoarthritis Uncover Underlying Heterogeneity Associated With Disease Progression

Introduction

Osteoarthritis (OA) represents the most frequent age-associated 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 (Li and Zheng, 2020). The immune system that is involved in the initiation and progression of OA is a key element in the pathogenesis of this disease (Zhang et al., 2019). The inflammatory microenvironment of OA synovial tissues consists of immune cells and inflammatory mediators such as inflammatory cytokines (IL-1β, TNFα, IL-6, IL-15, IL-17, and IL-18) 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.

Materials and Methods

OA Expression Profiles

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.

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

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

Western Blotting

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.

Transfection

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.

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-2′-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).

Results

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; 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).

FIGURE 1. Identification of DEGs for OA synovial tissues. (A) Volcano diagram of DEGs in OA compared to healthy synovial tissues in the GSE55235 dataset. Red: up-regulation and green: down-regulation. X-axis: −log10(adjusted p-value) and y-axis: log[fold-change (FC)]. (B) Heatmap of DEGs in OA (red) than healthy (blue) synovial samples in the GSE55235 dataset. N: normal and O: osteoarthritis. (C,D) Volcano and heatmap plots of DEGs in OA than healthy synovial tissue specimens in the GSE55457 dataset. (E) Venn diagram of common DEGs in the GSE55235 and GSE55457 datasets.

TABLE 1. 57 common DEGs for OA in the GSE55235 and GSE55457 datasets.

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

FIGURE 2. GO and pathway enrichment analysis of DEGs. (A) The top ten biological processes of common DEGs. X-axis represents the number of enriched DEGs and y-axis represents the biological process terms. The more the color tends to red, the smaller the adjusted p-value. (B) The top 30 KEGG pathways enriched by common DEGs.

TABLE 2. The detailed information of GO and pathway enrichment analysis results.

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

FIGURE 3. Landscape in immune infiltration between OA and healthy synovial tissue specimens from the GSE55235 and GSE55457 datasets. (A) PCA of synovial tissues based on immune cell infiltration levels. The first two principal components which decipher the most of the data variation. (B) The relative percent of 22 kinds of immune cells in each specimen. (C) Heatmap of the infiltration levels of immune cells in OA (red) and normal (blue) synovial tissues. (D) Positive (red) and negative correlations between different immune cells in synovial specimens. Correlation coefficient is marked in each box.

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

FIGURE 4. Comparisons of immune infiltration between OA and healthy synovial tissues. (A) Violin plots of the infiltration levels of immune cells in OA (red) and healthy (blue) synovial tissue specimens. Scatter plots for the significant differences in infiltration levels of (B) dendritic cells activated, (C) T cells CD4 memory resting, (D) neutrophils, (E) mast cells activated, (F) macrophages M0, (G) mast cells resting and (H) monocytes in OA and healthy synovial tissues.

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

FIGURE 5. Quality control of scRNA-seq of osteoarthritic synoviocytes in the GSE152805 dataset. (A) Violin plots of the number of feature genes, the count of total genes and the percent of mitochondrial genes after filtration. (B) Scatter plots of correlations of the count of total genes with the percent of mitochondrial genes and the number of feature genes. (C) The top 2000 highly variable genes and visualization of the top ten genes. (D,E) Visualization of dimensionality reduction results of the scaling data. (F) Elbow plot for identifying the optimal number of PCs.

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

FIGURE 6. Identification of cell clusters and marker genes in OA synoviocytes. (A) Visualization of cell clusters marked by different colors for OA synoviocytes by t-SNE method. (B) Heat map of the scaled expressions of the top ten marker-genes for each cell cluster. (C) Annotation of t-SNE cell types for OA synoviocytes. (D) Cell cycle (G1, G2M, S) of OA synoviocytes on the t-SNE map. (E–G) Visualization of the expression levels of cytokine receptor marker genes for each cluster. (H–J) Visualization of the expression levels of inflammatory cytokine marker genes for each cluster.

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

FIGURE 7. Identification of five immune cell subpopulations and associated marker genes for OA synoviocytes. (A) Visualization of five cell subpopulations with unique colors. (B) Heatmap for the scaled expression of the top 20 marker genes for each cell subpopulation. (C,D) Visualization of the expression levels of cytokine receptor marker genes for each cell subpopulation. (E,F) Visualization of the expression levels of inflammatory cytokine marker genes for each subpopulation. (G) Visualization of the expression levels of surface markers for each subpopulation.

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.

FIGURE 8. Pseudo-time analysis of osteoarthritic synoviocytes by Monocle 2 algorithm. (A) Distribution of osteoarthritic synoviocytes in 7 different states. (B) The mRNA expression of inflammatory cytokines and cytokine receptor in each state. (C) Distribution of osteoarthritic synoviocytes according to pseudo-time value. The darker the color, the larger the pseudo-time value. (D) Genes that dynamically change with pseudo-time value.

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.

FIGURE 9. Ligand-receptor network based on B cells, dendritic cells, macrophages, monocytes, and T cells.

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

FIGURE 10. Biological functions of differentially expressed marker genes in five immune cells. (A,B) Before and after batch effects of the integrated GSE55235 and GSE55457 datasets. (C) Volcano plots of differentially expressed genes between OA and normal synovial tissues. (D–J) Functional annotation analysis of differentially expressed marker genes in (D) B cells, (E,F) dendritic cells, (G) macrophages, (H,I) monocytes, and (J) T cells.

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

FIGURE 11. Silencing NR4A1 and NR4A2 relieves IL-1β-induced chondrocyte damage. (A–C) Western blotting of the expressions of NR4A1 and NR4A2 in IL-1β-induced chondrocytes. (D–F) Western blotting of the expressions of NR4A1 and NR4A2 in IL-1β-induced chondrocytes after silencing NR4A1 or NR4A2. (G–I) Flow cytometry detecting the apoptosis of IL-1β-induced chondrocytes with knockdown of NR4A1 and NR4A2. (J–L) EdU assay measuring the proliferation of IL-1β-induced chondrocytes with knockdown of NR4A1 and NR4A2. (M–T) Western blotting of the expressions of MMP3, MMP3 and Collagen II in IL-1β-induced chondrocytes under knockdown of NR4A1 and NR4A2. **p < 0.01; ***p < 0.001; ****p < 0.0001.

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, OA-related DEGs were also distinctly enriched in TNF pathway (Li et al., 2019). 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, Toll-like 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 (Li et al., 2020). 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.

Funding

This work was funded by Sun Yat-sen Clinical Research Cultivation Program of Sun Yat-sen Memorial Hospital, Sun Yat-sen University (SYS-Q-202105); Medical Research Foundation of Guangdong Province (A2021280).

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/fmolb.2021.748360/full#supplementary-material.

Supplementary Table S1 | DEGs of OA synovial tissues in the GSE55235 dataset.

Supplementary Table S2 | DEGs of OA synovial tissues in the GSE55457 dataset.

Abbreviations

DEGs, differently expressed genes; FC, fold-change; FDR, false discovery rate; GEO, gene expression omnibus; GO, gene ontology; KEGG, kyoto encyclopedia of genes and genomes; limma, linear models for microarray data; OA, osteoarthritis; PCA, principal component analysis; PCs, principal components; scRNA-seq, single-cell RNA sequencing; Tregs, T cells regulatory; t-SNE, t-distributed statistical neighbor embedding.

References