Single-cell RNA seq analysis identies the biomarkers and differentiation of chondrocyte in human osteoarthritis

Background (cid:0) Single-cell RNA sequencing (scRNA-seq) was recently adopted for exploring molecular programmes and lineage progression patterns of pathogenesis of important diseases. In this study, scRNA-seq was used to identify potential markers for chondrocytes in osteoarthritis (OA) and to explore the function of different types of chondrocytes in OA. Methods (cid:0) Here we aimed to identify the biomarkers and differentiation of chondrocyte by Single-cell RNA seq analysis. GeneOntology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were used to identify the function of candidate marker genes in chondrocytes. Protein–protein interaction (PPI) network was constructed to nd the hub genes in 3 types of chondrocyte respectively. We also used qRT-PCR to detect the expression level of the candidate marker genes in different types of chondrocyte. Results (cid:0) In this study, we characterized the single-cell expression proling of 480 chondrocyte samples and found hypertrophic chondrocyte (HTC), homeostatic chondrocyte (HomC) and brocartilage chondrocyte (FC) respectively. The results of GO and KEGG analysis showed the candidate marker genes made specic function in these chondrocytes to regulate the development of OAs respectively. We further revealed the differential expression of top 10 marker genes in 3 types of chondrocyte. The marker genes of HTC and FC were mainly expressed in their cell subset respectively. The marker genes of HomC did not have obviously differential expression among different types of chondrocyte. Last, we predicted the key genes in each cell subset. CD44, JUN and FN1 were predicted tightly related to the proliferation and differentiation of chondrocytes in OAs and could be regarded as biomarkers to estimate the development of OA. Conclusion (cid:0) Our results provide new insights into exploring the roles of different types of chondrocyte in OA. The biomarkers of chondrocyte were also valuable for estimating OA progression.

chondrocyte (FC) respectively. The results of GO and KEGG analysis showed the candidate marker genes made speci c function in these chondrocytes to regulate the development of OAs respectively. We further revealed the differential expression of top 10 marker genes in 3 types of chondrocyte. The marker genes of HTC and FC were mainly expressed in their cell subset respectively. The marker genes of HomC did not have obviously differential expression among different types of chondrocyte. Last, we predicted the key genes in each cell subset. CD44, JUN and FN1 were predicted tightly related to the proliferation and Background Osteoarthritis (OA) is the most common type of arthritis. It is a disability, delayed and degenerative disease characterized by appearance of in ammation of articular cartilage and synovium, resulting in stiffness, swelling, pain and loss of mobility [1,2]. An important pathological feature of OA is the degradation of articular cartilage. Cartilage is composed of chondrocytes, stroma and bers [3]. Chondrocyte is the only cell type present in mature cartilage and is responsible for repairing damaged cartilage tissue [4]. Dysfunction of chondrocytes will result in degradation of the matrix, cell death and total loss of cartilage integrity [5]. Some reports found chondrocytes in aging or diseased cartilage became senescent with associated phenotypic changes contributing to the development or progression of OA [6]. Apoptotic chondrocyte death also occurs more frequently in OA compared to normal cartilage and in OA lesional compared to OA non-lesional cartilage [7]. Owing to the important roles of chondrocytes in the pathogenesis of OA, some researchers have tried to determine the development of OA based on the differentiation and growth of chondrocytes [8,9]. A few biochemical markers have been found and were used for estimating growth stage of chondrocytes in OA. However, biochemical markers that can effectively predict cell subpopulations of chondrocytes remain largely unknown.
The formation of chondrocyte begins with the aggregation of mesenchymal cells [10]. Mesenchymal cells undergo a series of changes such as aggregation, migration, proliferation and fusion, and eventually differentiate into chondrocytes [11]. Chondrocytes undergo in four stages: resting stage, proliferation stage, prehypertrophy stage, and hypertrophy stage, and nally calcify into bone [12]. There exist some changes in gene and protein expression, cell morphology, cell surface markers and metabolic activity in the process of chondrocyte differentiation. Of which, the expression of some genes and proteins can not only regulate chondrogenic differentiation but also determine the stage of chondrogenic differentiation. Sox5 (L-Sox5), Sox6 and Sox9 have been used as biochemical markers of proliferative chondrocyte (proC) and cooperatively activated expression of the chondrocyte differentiation marker COL2A1 [13].
COL10A1 and ALP are markers of hypertrophic chondrocyte (HTC), when chondrocytes are transferred from proliferation stage to hypertrophy stage, the expression levels of these two genes are signi cantly increased [14]. Runx2 is mainly expressed in the perichondrium and prehypertrophic chondrocytes (preHTCs), which regulate the transformation of chondrocytes into HTCs [15]. These biomarkers play important roles in de ning types of different chondrocytes. However, because of the limited number of available markers to identify the cell subtype of chondrocytes in OA, the exact process of chondrocyte differentiation in OA is not clear.
In this study, we derived and characterized the cell features and marker genes among different cell subtype of chondrocytes using scRNA-seq pro ling from chondrocytes isolated from 3 patients with OA.
We conducted computational analysis to de ne the cell subtype of chondrocytes and to identify the Table1. We ltered out poor-quality cells at rst. The transcriptome data were merged into one matrix and the normalization process were conducted by using the R package of limma.
Processing of single-cell RNA-seq data We utilized R packages of limma, Seurat, dplyr and magrittr to lter out the poor data of matrix to transcriptome and normalize the data. We rst read the data of matrix and merged the repeated genes. We used R package of Seurat to generate the object and ltered out cells with poor quality [17]. The percentage of the gene numbers, cell counts and mitochondria sequencing count were calculated. Per Seurat's recommendation for quality control, the gene expressed less than 3 cells and the cell contained less than 50 genes were discarded. To remove potential cell doublets and low quality cells, the cells that have the percentage of mitochondrial genes more than 5% were ltered out [18]. We also extracted genes that have a higher coe cient of variation between cells and constructed a characteristic variance diagram.
After the quality control of data for the 3 OA patients and found the gene symbols, we performed a principal component analysis (PCA) analysis and the top principal component for the cells was obtained [19]. We estimated the top 20 principal component by P value and the signi cant components were selected for t-Distributed Stochastic Neighbor Embedding (t-SNE) algorithm. For visualization, we performed the t-SNE analysis based on that selected PC dimensions with the R package of Seurat [20]. For classifying all ltered cells, we set the clustering parameter resolution to 0.5 for the function FindClusters in Seurat. The standard of expression feature of marker genes was cutting off cell samples and marker genes with log fold change lter=0.5 and adj PvalFilter=0.05 by the R package of limma [21]. The R package of ggplot2 was used to nd the candidate marker genes in different cluster. The clusters and cell categories were also annotated by candidate marker genes on the R package of SingleR. Monocle pseudotime trajectory of cells were also constructed by the R package of monocle [22].
Culturing the ATDC5 cells and L-929 cells ATDC5 cells were used as a vitro model for chondrocytes. It is a mouse chondrogenic cell line which showed a sequential transition of phenotype in vitro. ATDC5 has been derived from teratocarcinoma AT805. ATDC5 is an important chondrocyte model, which has been widely applied to the study of chondrocyte proliferation and differentiation. The ATDC5 cells (1×10 4 cells, purchased from ShZeYe, Shanghai, China) were inoculated in a cell vial containing DMEM/F12 (1:1) and 5% FBS. When the cells were more than half of the wall, added 1% of ITS (Insulin, Transferrin, Selenium). Changed the cell culture medium every 2d. The chondrocytes of 7d,14d and 21d were selected to identify the growth stage. The chondrocytes were divided into proC, HTC and resting chondrocyte (resC) by the expression of marker genes for different phase. All cells were cultured in a constant temperature incubator at 37°C with 5% CO 2 .
RNA extraction and reverse-transcription quantitative PCR (RT-qPCR) Total RNA from tissues was extracted and was reverse-transcribed into cDNA using RNeasy Mini kit (Qiagen, Valencia, CA, USA) according to the manufacturer's instructions. The expression of mRNAs was determined using SYBR Premix Ex Taq TM II (Takara Bio Inc., Japan). The primers using GAPDH as internal parameter. The relative expression level was calculated using the 2 -ΔΔCt method. The experiment was repeated 3 times. The primers were listed in Supplementary Table 2 Protein-protein interaction (PPI) network analysis Gene expression data from 3 clusters were rst mapped with Retrieval of Interacting Genes (STRING) database (http://www.string-db.org/) and then constructed PPI network. Then, the interaction data were then put into Cytoscape software to analyze the degree of connectivity. The top 20 candidate marker genes with high degree of connectivity were selected as the hub genes and visualized in histogram.

Statistical analysis
Statistical analysis was performed using R packages. Bar graphs show average expression levels ± SEMs. P values were calculated using a two-tailed Student's t-test. All statistical calculations were processed using SPSS 22.0 statistical software. The evaluation of RT-qPCR data was performed using one-way ANOVA with Tukey's post hoc test. P < 0.05 was considered signi cant.

Single-cell pro ling of human OA cartilage chondrocytes
To explore the genes closely related to the development of chondrocytes, we acquired 480 chondrocyte samples from 3 patients with OA (OA1, OA2 and OA3). Each OA included 160 chondrocytes. The program Seurat were used to analyze data. The distribution pattern chart showed the gene numbers and the sequencing count of each cell ( Figure 1a). Most cells of the 3OAs included 1000-5000 genes. The sequencing count of cells focus on the range of 950000-980000. We did not nd the mitochondrial sequencing count from the data. By analyzed the correlation between nFeature and nCount, the signi cantly positive correlation between nFeature and nCount were found ( Figure 1b). We also identi ed the gene symbols with signi cant difference across cells (Figure 1c), of which10 signi cantly differentially expressed genes (COL1A1, ZNF585B, SLC5A12, COL1A2, IBSP, SPARCL1, TMSB4X, SPP1, CCL2, COL14A1) across the cell samples have been found. with an estimated P value, and we selected the signi cant components for subsequent analysis ( Figure  1e). Subsequently, cells were clustered based on a graph-based clustering approach, and were visualized in dimensional space using t-distributed Stochastic Neighbor Embedding (tSNE). By tSNE algorithm we classi ed the cells into 6 clusters ( Figure 1f). However, subgroups of the cells could not be distinguished.
The 6 clusters were all annotated as chondrocyte. Therefore, a deeper analysis was needed to nd the subgroups of cells.

Cell trajectory analysis and annotations
In order to nd the subgroups of cells, we rst revealed the tendency curve of different clusters by cell trajectory analysis. Before conducting the cell trajectory analysis, we performed differential analysis using the limma package and identi ed a total of 1010 candidate marker genes with | log fold change | > 0.5 and adj Pval < 0.05 (Supplementary Table 3). The top 10 differential genes between the six clusters in the heatmap plot were illustrated in Figure 2a. Moreover, we constructed the trajectory analysis to reveal the tendency curve of the six clusters (Figure 2b). The results showed the cells of cluster2, cluster3 and cluster5 may had signi cant difference between each other and may belonged to different subgroup of chondrocyte. In addition, we empirically de ned the subgroup of chondrocyte in each cluster by candidate marker genes. By annotation, we found the cells of cluster2 belonged to HTCs. The cells of cluster3 belonged to homeostatic chondrocytes (HomCs) which were de ned by QB and his partners [16]. The cells of cluster5 belonged to FCs. We failed to annotate the cells of cluster0, cluster1 and cluster4. We suspected that the cells of cluster0, cluster1 and cluster4 may belong to different period of HTCs, HomCs and FCs. However, these existed a low degree of heterogeneity of the candidate marker genes among the 3 clusters. Therefore, we could not de ne the speci c cell subtype or cell period of the 3 clusters. By analyzing these results, we suggested the trajectory analysis combined with the annotation to candidate marker genes of tSNE was helpful to de ning the types of cells between different clusters. For the top 10 candidate marker genes of HTC (cluster2), HomC (cluster3) and FC (cluster5), we illustrated the signi cantly differential expressions of these genes in six clusters via bubble plot (Figure 2c-e). The plots showed most of these genes were speci cally overexpressed in the cells of each cluster.

Function and pathway enrichment analysis
By cell trajectory analysis and annotation of candidate marker genes in 3 cluster (cluster2, cluster3 and cluster5), we de ned HTCs, HomCs and FCs. To testify the validity of the cell annotation to 3 cluster and the function of candidate marker genes in the cells of 3 cluster, we carried out the GeneOntology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of the candidate marker genes in the 3 cluster respectively. For the genes belonging to HTCs, most functions including oxidative stress, epithelial cell proliferation and growth factor activity were closely related to HTCs (Figure 3a). Similar to GO analysis, the KEGG analysis of the genes in HTCs also had some important pathways (ECM-receptor interaction, IL-17, TNF) ( Figure 3b). The genes of these pathways all play important roles in HTCs. The results of GO and KEGG analysis to genes in HomCs showed some functions such as unfolded protein, temperature stimulus and heat shock protein could modulate cellular homeostasis (Figure 3c-d). These functions and pathways illustrated HomCs may regulated cell cycle, development, RNA metabolism and biosynthesis. FCs have close relationship with collagen formation. The function and pathway enrichment analysis of the genes in cluster 5 showed some notably functions with collagen (collagen metabolic, collagen trimer and collagen binding) (Figure 3e-f). These results further evidence the correctness of the annotation to cells in 3 clusters. The results of GO and KEGG analysis also indicate these candidate marker genes may play signi cant roles during the process of proliferation and differentiation for chondrocytes.
Determining the expression level of candidate marker genes among resCs, proCs, HTCs and FCs The GO and KEGG analysis identi ed the annotation of clusters was accurate. In order to con rm the candidate marker genes indeed had aberrant expression in different subtypes of chondrocyte, we detected the expression level of top 10 candidate marker genes in 3 clusters among resCs, proCs, HTCs and FCs. We rst used ATDC5 cells to separate resCs, proCs and HTCs. We separated the 3 subtypes of chondrocyte by comparing the differential expression of marker genes (resC: SOX9 and Fn1, proC: COL2A1 and Runx2, HTC: COL10A1 and COL1A1) among 7d, 14d and 21d in ATDC5 cells (Supplementary Figure 2). For FCs, we used L-929 cells which were proximal cartilage broblasts of mouse as FCs.
After nishing the preparation of subtypes of chondrocyte, we detected the expression level of top 10 candidate marker genes in 3 clusters among resCs, proCs, HTCs and L-929 cells. For HTCs, the results showed the 10 genes all had higher expression level in recCs, proC, HTCs compared with L-929 cells. The genes of TNFAIP6, GAS1, SLC7A2 and TNC were markedly upregulated in HTCs (Figure4a-b). The results of top 10 genes in HomCs showed the 10 genes all had higher expression level in recCs, proC, HTCs. However, compared with top 10 genes in HTCs the expression level of most genes in HomCs did not have obviously differential expression among resCs, proCs, HTCs (Figure4c-d). When examined the expression level of top 10 genes in FCs, we found the expression of most genes up-regulated in L-929 except for CFD, SFRP2 and SPARCL1. The gene expression of COL1A1 and COL1A2 also showed up-regulated in HTC compared with resC and proC (Figure4e-f). These results implied that most candidate marker genes may have special expression in their own subgroup.
Identi cation of hub genes from candidate marker genes in HTCs, HomCs and FCs To further identify the key genes in each cell subset. Protein-protein interaction (PPI) network of candidate marker genes in 3 clusters was constructed by Cytoscape software respectively, based on STRING database (Supplementary Figure 3). The number of genes which were used to construct PPI in HTCs, HomCs and FCs were 104, 175 and 329, respectively. The top 20 candidate marker genes (HTCs: the degree of connectivity≥8, HomCs: the degree of connectivity≥19, FCs: the degree of connectivity≥40) with high degree of connectivity were displayed (Figure5a-c). These top genes might play critical roles in chondrocyte progression.
CD44, JUN and FN1 were the hub genes with highest degree of connectivity in each cell subset. In order to test the relationship between candidate marker genes and cell subset. The expression of the 3 genes in the 4 chondrocyte subsets was detected. The expression of CD44 was gradually increased from resC to HTC of chondrocytes (Figure 5d). The result of bubble plot also showed CD44 was speci cally overexpressed in HTCs (cluster2) (Supplementary Figure 4). JUN was a candidate marker gene that was speci cally overexpressed in HomCs (cluster3). The expression of JUN was contrary to CD44, it was high expressed in resC, proC and HTC compared with L-929 (Figure 5d). The expression level of FN1 was markedly upregulated in L-929 which was consistent with the result of FN1 in bubble plot (Supplementary Figure 4). These results combined with previous analysis con rmed some candidate marker genes can reveal the progression of chondrocytes. CD44, JUN and FN1 may play important roles in different subtypes of chondrocyte.

Discussion
OA is the most common disease of joints and currently lacks treatment options that modify structural pathology [23,24]. It is usually de ned as the result of an imbalance in the degradation and synthesis of the normal coupling of chondrocytes, cell matrix, and subchondral bone under the in uence of mechanical and biological factors [24,25]. It is now generally accepted that the chondrocytes are the target of these abnormal biomechanical factors, and that biochemical and genetic factors also contribute to alterations in the normal functional activities of these cells [26]. Clinically, some researchers can judge the progress of osteoarthritis patients by the degradation of cartilage tissue proteoglycan and collagen bers [27,28]. Several genes have also been found to play important roles in the development of chondrogenesis in osteoarthritis and can serve as markers for the development of osteoarthritis. JK and his partners found lncRNA FAS-AS1 could promote the degradation of extracellular matrix of cartilage in osteoarthritis by increasing the expressions of MMP1 and MMP13 and decreasing the expression of COL2A1 in chondrocytes [28]. These genes could also be used for predicting the internal state of chondrocytes in OA pathogenesis. However, due to the complexity of chondrocyte differentiation in osteoarthritis and lack of speci c biomarker to characterize chondrocyte differentiation in osteoarthritis, the differentiation process and cell subtype of chondrocytes in osteoarthritis remain largely unknown.
Here, we investigated OA chondrocytes at different stages by using scRNA-seq analysis. We found 3 populations of chondrocytes (HTCs, HomCs and FCs), new markers of chondrocyte populations and signaling pathways involved in OA pathogenesis based on scRNA-seq analysis. We also found the hub genes which may play important roles in the 3 chondrocyte populations respectively. There results may contribute to the early diagnosis and therapeutic treatment of human OA.
For chondrocyte taxonomy, some studies have shown chondrocyte subtypes mainly included ProCs, preHTCs and HTCs [29][30][31]. In addition to these three common chondrocyte subtypes, QB and his partners rst identi ed 4 new populations of chondrocytes in the human OA cartilage by performing unbiased transcriptome-wide scRNA-seq analysis and computational analysis on some chondrocytes from 10 OA patients [16]. QB and his partners constructed a single-cell analysis method to separate and identify chondrocyte populations at single-cell resolution and at the whole-transcriptome scale. The datasets and results of their report could help researchers to de ne the function of chondrocyte subtypes and nd new for chondrocyte subtypes. However, the analytical method of their study was extremely complex. For our research, low-Volume cell samples were used to de ne the chondrocyte subtypes and identify novel markers of chondrocyte subtypes by general analytical method for scRNA-seq. In this study, we found 3 cell populations (HTCs, HomCs and FCs) in chondrocytes of 3 OAs. HTCs are divided into prohypertrophy and hypertrophy [32]. In recent years, the research on HTCs has become a hot topic in the research of osteoarthritis. Many HTCs were found in the degenerative articular cartilage [33]. These chondrocytes have different types of collagen which have different roles in the development of osteoarthritis. It is important to study the speci c role of HTCs in osteoarthritis and the genes that regulate chondrocytes. In our study, we also found HTCs existed in the chondrocytes of OAs. The GO analysis of HTCs concentrated on oxidative stress, epithelial cell proliferation and growth factor. Some researchers found SIN-1 could induce oxidative stress lead to necrosis in HTCs [34,31]. For cell proliferation and growth factor, some growth factor such as COL10A1 and ALP were high expression in HTCs and promoted the proliferation of HTCs [35,36]. Therefore, the candidate marker genes of HTCs may have close relationships with proliferation and differentiation of HTCs in OAs. The expression analysis of top 10 candidate marker genes in different types of chondrocyte also showed TNFAIP6, GAS1, SLC7A2 and TNC were high expression in HTCs. TNFAIP6, GAS1 and TNC have been found associated with the biological function of chondrocyte or cartilage by other researchers [37][38][39]. The change in expression level of these genes in HTCs of OAs can be used to judge the development status of OA.
Compared with HTCs, there are few reports on the role of HomCs and FCs in osteoarthritis. FCs are new cell types that have been discovered to differentiate from other chondrocyte subset in recent years. FCs was mainly found in the late stages of OA and possessed a high ratio of genes related to unfavorable OA outcomes and a capacity for vascularization [16,40]. The GO and KEGG analysis showed FCs have close relationship with collagen formation. The macromolecular framework of articular cartilage consists of collagens (predominantly type II collagen), proteoglycans, and non-collagenous proteins [3,41]. Of which, collagen is the matrix that provides the cartilage with its tensile stiffness and strength. It is possible for researchers to nd new treatments of OA by regulating the collagen formation and proliferation of FCs which could secret collagen. We also found 6 candidate marker genes (S100A4, COL1A2, TMSB4X, THY1, COL14A1 and IFI27) were high expression in FCs. These genes may be used as biomarkers for distinguishing FCs from other chondrocytes. HomCs were helpful to regulating the homeostatic balance of chondrocyte metabolic activity in response to environmental signals. This type of chondrocyte may play an important role in maintaining rhythmic behaviors and daily cycles of metabolism during OA cartilage degeneration [42,43].Our research also found HomCs may played important roles in regulating cell cycle, RNA metabolism and biosynthesis. We also found all top 10 candidate marker genes of HomCs were high expression among resCs, proCs, HTCs. These results imply these candidate marker genes of HomCs might regulate the circadian clock system of chondrocyte in OAs.
Three hub genes were also found among these candidate genes in HTCs, HomCs and FCs respectively. CD44 is the membrane receptor for hyaluronan. CD44 is the primary receptor of retention and anchoring of proteoglycan aggregates to the chondrocyte plasma membrane [44,45]. This receptor also play important roles in cartilage such as maintains cartilage homeostasis and affects both chondrocyte survival pathways as well as apoptotic pathways [46,47]. In this study we regarded CD44 as the hub gene of HTCs and found this receptor was high expression in HTCs. These results show CD44 is a biomarker of HTCs in OA and this receptor exerts more in uence on HTCs than other types of chondrocyte in OA. JUN is also known as AP-1. The abnormal expression of JUN in chondrocyte can affect the differentiation of chondrocyte in growth cartilage. MH and his partners found speci c inhibition of c-Fos/AP-1 could regulate differentiation of chondrocyte in OA which effectively prevented cartilage destruction and osteophyte formation [48]. Some researchers also found the suppression of c-Fos/AP-1 expression could regulate the activity of OA chondrocyte under pathological conditions [49]. In our study, we found JUN which was the hub gene of HomCs was high expression among different types of chondrocyte except for FCs. The function of HomCs mainly focus on maintain rhythmic behaviors and daily cycles of metabolism during OA cartilage degeneration. We speculate JUN may inhibit the differentiation of chondrocyte by affecting cell cycle of chondrocyte. Fibronectin1 is a glycoprotein present in a soluble dimeric form in plasma, and in a dimeric or multimeric form at the cell surface and in extracellular matrix. Some researchers found FN1 could increase chondrocyte cell death [50]. The high expression of FNI in FCs in our study implied this gene could be used for testing the growth of brocartilage.
Nevertheless, there were still several weaknesses for further optimization. First, the method of scRNA-seq analysis in our research may only applied to low-Volume cell samples. Compared with the analysis method in the study of QB and his partners, we also need to optimize the scRNA-seq analysis method to focus on large-Volume Sample stacking [16]. Second, the trajectory analysis combined with empirically annotation of the chondrocyte subsets could de ne the speci c cell subtype of 3 clusters. However, cluster0, cluster1 and cluster4 were still unde ned. The 3 clusters should be performed a deeper annotation with more reasonable method in the future.

Conclusions
In conclusion, our single-cell analysis identi ed 3 types of chondrocyte with their candidate marker genes in OA. Our analysis identi ed the function and discriminative markers of HTCs, HomCs and FCs in OA.
The expression level of hub genes and the growth of 3 types of chondrocyte have been predicted with close relationships with the development of OA. These biomarkers may provide important diagnostic and therapeutic strategies for patients with OA in the future.
Declarations Acknowledgments Not applicable.
Authors' contributions XL, NL, and RF were involved in acquisition and interpretation of the data. XL and SQ supervised cohort construction and contributed to the design of the work. LM and QF analyzed and interpreted data. LQ and XL cultured cells and conducted the experiments of RT-qPCR. XZ and SQ wrote the paper. XL revised the text and is the guarantor of this work. All the authors have read and approved the manuscript. All authors have approved the submitted version and have agreed both to be personally accountable for the author's own contributions and to ensure that questions related to the accuracy or integrity of any part of the work, even ones in which the author was not personally involved, are appropriately investigated, resolved, and the resolution documented in the literature.

Funding
The study was funded by Natural Science Foundation of Fujian Province in China through the project "The inhibition of TNF-α on osteolysis around the prosthesis after joint replacement in patients with rheumatoid arthritis" to XL (2019J01050826) Availability of data and materials The  we illustrated the signi cantly differential expressions of FN1, JUN and CD44 in six clusters via bubble plot.
Supplementary Table 1 Clinical and demographic characteristics of the study population.
Supplementary   Cell trajectory analysis and annotation of 6 clusters a: The heatmap of top 10 cluster speci c genes from each cluster. Yellow -high expression. Purple -low/no expression. Each row represents a gene, and each column represents a single cell; b: Cell trajectory analysis revealing the OA chondrocyte progression colored according to clusters. c-e: The expression levels of 10 candidate marker genes in HTC (c), HomC (d) and FC (e) among six clusters were illustrated via bubble plot.

Figure 3
Function and pathway enrichment analysis of the candidate marker genes in HTCs, HomCs and FCs. a-b: GO and KEGG analysis of the candidate marker genes in HTCs; c-d: GO and KEGG analysis of the candidate marker genes in HomCs; e-f: GO and KEGG analysis of the candidate marker genes in FCs;

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.