Comparison of the major cell populations among osteoarthritis, Kashin–Beck disease and healthy chondrocytes by single-cell RNA-seq analysis

Chondrocytes are the key target cells of the cartilage degeneration that occurs in Kashin–Beck disease (KBD) and osteoarthritis (OA). However, the heterogeneity of articular cartilage cell types present in KBD and OA patients and healthy controls is still unknown, which has prevented the study of the pathophysiology of the mechanisms underlying the roles of different populations of chondrocytes in the processes leading to KBD and OA. Here, we aimed to identify the transcriptional programmes and all major cell populations in patients with KBD, patients with OA and healthy controls to identify the markers that discriminate among chondrocytes in these three groups. Single-cell RNA sequencing was performed to identify chondrocyte populations and their gene signatures in KBD, OA and healthy cells to investigate their differences as related to the pathogenetic mechanisms of these two osteochondral diseases. We performed immunohistochemistry and quantitative reverse-transcription PCR (qRT-PCR) assays to validate the markers for chondrocyte population. Ten clusters were labelled by cell type according to the expression of previously described markers, and one novel population was identified according to the expression of a new set of markers. The homeostatic and mitochondrial chondrocyte populations, which were identified by the expression of the unknown markers MT1X and MT2A and MT-ND1 and MT-ATP6, were markedly expanded in KBD. The regulatory chondrocyte population, identified by the expression of CHI3L1, was markedly expanded in OA. Our study allows us to better understand the heterogeneity of chondrocytes in KBD and OA and provides new evidence of differences in the pathogenetic mechanisms between these two diseases.


Introduction
Kashin-Beck disease (KBD), an endemic and chronic degenerative osteochondral disease with irreversible pathological and clinical signs, is characterised with clinical manifestations including shortened and enlarged fingers at early stage of onset of KBD, and deformed limb joints, limited movement and even dwarfism in some patients with advanced disease 1,2 . Osteoarthritis (OA) is the most common chronic disease, involving progressive joint dysfunction with ageing. Several similar pathological changes have been found in KBD and OA, including extracellular matrix (ECM) degradation, cartilage lesions, and reduction and disruption of proteoglycans (PGs) and collagens [3][4][5][6] . However, the age of onset, symptoms and Xray findings differ between KBD and OA 3 , which suggests that they are caused by different aetiologies and have distinct pathogeneses. Meanwhile, adult KBD patients typically have advanced KBD, often accompanied by OA 4,6 . Recent studies have found differences in various cellular pathways between KBD and OA chondrocytes 5,7,8 , suggesting differences in multiple processes, such as metabolism, apoptosis, adaptive immune defence, cytoskeleton, cell movement and extracellular matrix turnover 4,[9][10][11][12] . However, these studies do not take into consideration the heterogeneity of articular cartilage cell types and therefore may be missing important information on cartilage degeneration processes specific to KBD or OA.
Articular cartilage is a highly organised structure composed of multiple zones, namely, the superficial, middle, deep and calcified zones, which display three chondrocyte phenotypes: persistent, transient and hypertrophic cells 13,14 . The chondrocyte phenotype, cell shape, and extracellular matrix (ECM) structure are different among the different zones [13][14][15] . In a previous study, seven defined populations of chondrocytes were identified in OA cartilage, specifically, proliferative chondrocytes (ProCs), prehypertrophic chondrocytes (PreHTCs), hypertrophic chondrocytes (HTCs), fibrocartilage chondrocytes (FCs) and three novel populations with distinct functions 16 . In addition, senescent cells (SNCs) and cartilage progenitor cells (CPCs) were also identified in recent studies [17][18][19][20] . Specific chondrocyte populations may play important roles in osteochondral diseases such as KBD and OA. For example, PreHTCs are capable of regulating the onset of hypertrophic differentiation, while HTCs can modulate the mineralisation of the surrounding matrix in the cartilage 16 , processes that are known to be more active in OA cartilage 21,22 . Additionally, CPCs have the ability to selfrenew and differentiate along multiple lineages, thereby contributing to cartilage repair and homeostasis, and their dysregulation may affect and modulate cartilage loss in the processes of OA and KBD 17,20 .
One of the main differences in pathological changes between KBD and OA is in the location of the primary cartilage loss. OA is characterised primarily by progressive degradation starting in the superficial zone, with subsequent inflammation. In contrast, the primary characteristic pathological changes of KBD patients were chondrocyte necrosis occurring in the deep zone of articular cartilage and hypertrophic layers of epiphyseal plate cartilage 23,24 . However, the pathological changes in the different chondrocyte populations in KBD cartilage have not yet been studied, and we therefore do not yet have a full understanding of the difference in cartilage degeneration mechanisms between OA and KBD.
Single-cell RNA sequencing (scRNA-seq), a novel and powerful method for investigating transcriptomic cell-tocell variation, is widely applied to identify various cell types and provide insights into the pathological processes of diseases [25][26][27][28] . In this study, scRNA-seq was used to better understand the changes occurring in the chondrocytes from patients with KBD and OA. The results allowed us to better understand the heterogeneity of chondrocytes between KBD and OA and provided us with new evidence to explore the differences in the pathogenetic mechanisms between these two diseases.

Single-cell trajectory of normal chondrocyte differentiation
To study the differentiation of healthy human chondrocytes into different clusters and the corresponding gene expression, we used monocle to reconstruct the pseudotime trajectory (Fig. 1B). We found that the root of the trajectory was mainly populated by CPCs-1 and CPCs-2, while FC-1 was distributed in the middle of the trajectory, FC-2 was distributed along the trajectory, RegCs were located in front of FCs, and PreHTCs were behind FCs. PreHTCs and HTCs were mainly distributed at the end of the trajectory. The two termini of the tree were populated by FC-1 and RegCs for fate 1 and FC-2, PreHTCs and HomCs for fate 2 (Fig. 1B).
Next, we assessed the expression of genes whose levels varied during CPC differentiation in cells at fates 1 and 2 of the trajectory. The expression of COL1A1, IGFBP6, CALR, SULF1, ACTB, LGALS1 and PFN1 was found to be similar in both groups. COL1A1, IGFBP6, CALR and SULF1 expression was slightly reduced, while ACTB, LGALS1 and PFN1 expression was slightly increased, at the early stage of differentiation, notably downregulated in fate 2 cells and upregulated in fate 1 cells. FMOD, DDIT4, IGFBP5, BGN, ZFAS1, FTL, COMP and CHI3L1 were slightly downregulated at the early stage of differentiation, notably upregulated in fate 2 and downregulated in fate 1. In contrast, CEMIP, ACTA2, PTX3, TAGLN, CCL2, ADAMTS1 and SERPINE1 expression decreased slightly in cells from the root to fate 2 but markedly increased in cells that differentiated via fate 1 ( Supplementary Fig. 1).
To study the distribution of different cell clusters, we used immunohistochemistry to detect marker gene expression (Fig. 1E). The results revealed that FCs were mainly distributed in the superficial and middle zones, CPCs were mainly distributed in the middle and deep zones, RegCs were mainly distributed in the superficial zone, and PreHTCs were mainly distributed in the superficial and deep zones.
Systemic comparison of the single-cell landscape among KBD, OA and healthy human chondrocytes We compared scRNA-seq among healthy chondrocytes, KBD chondrocytes and OA chondrocytes. In total, 16,375 cells were analysed, with 6140 cells from patients with KBD, 5834 cells from patients with OA and 4401 cells from healthy controls. Ten clusters were labelled by cell type using the expression of previously described markers, and one novel population, mitochondrial chondrocytes (MTCs), was identified by the expression of novel markers ( Fig. 2A). PreHTCs, identified as positive for S100A4, COL9A3 and AKR1C1, were divided into three primary subgroups of 6 clusters: the first expressing S100A4, TAGLN and SH3BGRL3 (S100A4 hi ), the second expressing COL9A3, SOX9 and COL11A1 (COL9A3 hi ) and the third expressing AKR1C1, LUM and AKR1C2 (AKR1C1 hi ). RegCs were found in six separate clusters, which were divided into two primary subgroups: the first expressing CTNNB1, TUBA1A and CDC42 (CTNNB1 hi ) and the second expressing CHI3L1, AEBP1 and STEAP1 (CHI3L1 hi ). MTCs, the novel population in this study, were found in two separate clusters, the first expressing MT-ND1, MT-ATP6 and MT-ND2 (MT-ND1 hi ) and the second expressing MT-CO1, MT-CO3 and MT-CYB (MT-CO1 hi ). FCs were found in seven separate clusters, which were divided into two primary subgroups: the first expressing IGFBP5, TPM1 and ADAMTS1 (IGFBP5 hi ) and the second expressing COL14A1, NOTCH3 and COL6A3 (COL14A1 hi ) (Fig. 2B, C and Supplementary  Fig. 2).
When separating t-SNE plots for each disease state, we noted that cluster 8, representing CPC cells, was present only in normal cartilage and was largely missing in both OA and KBD samples (Fig. 2D). In addition, an analysis of the proportion of total cells present in each population by sample and disease status revealed several significant changes among normal chondrocytes, KBD chondrocytes and OA chondrocytes (Fig. 2E). The populations of RegCs and prehypertrophic chondrocytes increased in OA chondrocytes compared with healthy control and KBD chondrocytes. Meanwhile, the numbers of hypertrophic chondrocytes, FCs and MTCs tended to increase from controls and OA to KBD.

Expanded cell population cluster in KBD samples
Two distinct chondrocyte populations seem to be expanded mainly in KBD samples compared to normal and OA samples. These include the novel MTCs and HomCs (Fig. 2D). HomCs were identified by their expression of the known markers MT1X, MT2A, MT1E, ATF5 and DDIT3. Another distinct population that was markedly expanded in KBD was the MTCs, which is a newly defined population that is identified by a signature that includes genes related to mitochondrial electron transport and the response to hydroperoxide, including the MT-ND1, MT-ATP6, MT-ND2 and MT-ND3 genes. The expression of MT-ND1 and MT1X in cartilage tissues was investigated by using immunohistochemistry assays, which validated the existence of MTCs and HomCs in KBD cartilage. The results showed that MTCs were mainly distributed in the middle and deep zones, and HomCs were mainly distributed in the superficial zone in KBD cartilage (Fig. 2F). Moreover, we validated the gene expression of MT-ND1 and MT1X in patients with KBD using RT-PCR and found that the results were consistent with the single-cell RNA-seq analysis ( Supplementary  Fig. 3).

Expanded cell population cluster in OA and normal samples
A distinct regulatory chondrocyte population that was markedly expanded in OA ( Fig. 2D) was identified by its expression of the known markers CHI3L1, AEBP1, PLIN2 and STEAP1 (Fig. 3A). The expression of AEBP1, CHI3L1 and STEAP1 in cartilage tissues was investigated by using immunohistochemistry assays and validated the existence of RegCs in OA cartilage. The results showed that AEBP1 was upregulated in the superficial and middle zones, and CHI3L1 and STEAP1 were upregulated in the superficial, middle and deep zones of OA cartilage tissues compared to normal cartilage tissues (Fig. 3B). In addition, we validated the gene expression of AEBP1 and CHI3L1 in patients with OA using RT-PCR and found that the results were consistent with the single-cell RNA-seq analysis ( Supplementary Fig. 3).
Distinct CPCs, which were markedly expanded in normal tissue (Fig. 2D), were identified by their expression of the known markers STMN1, PTTG1, CDKN3, BIRC5, UBE2S and CENPW (Fig. 4A). Although none of these markers were exclusive to CPCs, this cluster was the only population expressing all of these genes. We investigated the expression of CENPW and PTTG1 in cartilage tissues by using immunohistochemistry assays and validated the existence of CPCs in normal cartilage (Fig. 4B). The results revealed that CENPW was upregulated in the middle zone and PTTG1 was upregulated in the middle and deep zones in cartilage tissues from normal controls. In addition, the expression of PTTG1 and BIRC5 was investigated in cartilage tissues from OA and KBD patients, and the results showed that PPTG1 and BIRC5 were upregulated in the superficial zone of OA cartilage, while they were upregulated in the deep zone of KBD cartilage ( Supplementary Fig. 4).

Alignment of single-cell trajectory in normal, KBD and OA chondrocytes
To further investigate the relationship between the different disease statuses, we used monocle to study the transitions across specific cell types and sample states. We found that the pseudospace trajectory axis derived from monocle yielded three different trajectory states (Fig. 4C). Chondrocytes from normal cartilage were present at the start of the pseudospace trajectory, which then branched into 2 distinct fates, populated mainly with either KBD (fate 1) or OA (mostly fate 2) (Fig. 4D).
Next, we assessed the expression of genes regulated during chondrocyte differentiation in cells at fates 1 and 2 of the trajectory. The expression levels of FN1, CHI3L2, RGCC, COL9A3 and MDFI were markedly upregulated in cells at fate 1, whereas those of COL11A1, SMOC2 and FBXO2 were slightly downregulated at the early stage of differentiation and notably upregulated in cells at fate 1. HMOX1, COL1A1, WISP2, CD9 and PCOLCE were notably upregulated in cells at fate 2 and downregulated in cells at fate 1. IGF1, VCAN, PCSK5 and SPP1 were upregulated only in cells at fate 2 (Fig. 4E). These genes were differentially expressed at the prebranch point, suggesting that these genes may play a critical role in regulating the phenotype of chondrocytes at the beginning of the process of differentiation or development of chondrocytes from healthy controls to KBD or OA.

Fibrocartilage chondrocyte subpopulation
Upon examining FC populations in the combined analysis of KBD, OA and healthy human chondrocytes, we again identified distinct populations of IGFBP5 hi and COL14A1 hi FCs, each containing KBD, OA and healthy human chondrocytes. Two major and three minor subpopulations of FCs emerged (Fig. 6A, B). The first major population was defined by the expression of RPS4Y1 and PENK (RPS4Y1 hi FCs). A second major population was defined by the expression of COL3A1 (COL3A1 hi FCs). Three minor populations were distinguished by the expression of FGF1, IGFBP6 and COL14A1 (Fig. 6C).

Identification of effector and hypertrophic chondrocyte populations in KBD and OA chondrocytes
We separately analysed chondrocytes from patients with KBD and OA. Seven populations of chondrocytes were identified in both KBD and OA. The EC population was defined by the expression of MIA, S100B, FRZB, C2orf82 and S100A1 in both KBD and OA ( Supplementary Fig. 5), suggesting that these markers from ECs were characterised as having the same expression pattern in KBD and OA. In addition, HTCs were defined by the expression of IBSP, MMP13, RUNX2 and TMEM119 in KBD chondrocytes ( Supplementary Fig. 6).

Discussion
The pathogenesis of KBD and OA is different in many aspects. KBD mainly invades the epiphyseal cartilage, epiphyseal plate cartilage and articular cartilage and occurs in children aged 3-12 years. This progressive lesion seriously affects the growth and development of children's bones. Unlike KBD, OA is a degenerative chronic joint disease that mainly occurs in elderly individuals. Currently, we cannot distinguish KBD and OA well at the molecular level or propose precise treatment strategies for the two diseases. There are limited specific cell markers to report the internal status of chondrocytes in KBD and OA pathogenesis. Therefore, in this study, we performed a comprehensive analysis of  chondrocyte populations in KBD, OA and healthy control cartilage. We identified a new population of chondrocytes, MTCs, and, identified the relationships among KBD, OA and healthy control chondrocytes at the cellular level. The HomC and MTC populations were markedly expanded in KBD. The RegC population was markedly expanded in OA, and CPCs were markedly expanded in healthy controls. Our results thus identified the major cell populations in patients with KBD, OA and healthy controls and the corresponding marker genes, which provide a novel foundation for further exploration of the differences in the pathogenetic mechanisms involved in KBD and OA.
Currently, mitochondrial studies are involved in biomedicine due to their critical roles in ageing and disease development in humans 29 . Mitochondria are organelles that can generate most of the energy essential to the cell by converting nutritional molecules into adenosine triphosphate (ATP). An increase the production of reactive oxygen species (ROS) and a decrease in the consumption of ATP and oxygen by a series of metabolic alterations can be caused by mitochondrial dysfunction. Mitochondrial dysfunction also leads to an inflammatory response, which induces the synthesis of cytokines and matrix metalloproteinases (MMPs) 30 . Compared with that in normal cartilage, mitochondrial biogenesis in KBD and OA chondrocytes is altered, such that there is a reduction in the activities of complexes II, III, IV and V 31,32 , which causes failure in energy generation. Hence, mitochondrial function needs correction to prevent excessive production of ROS and oxidative stress 30 . In this study, we identified a new population of MTCs that are markedly expanded in KBD and express genes related to ATP synthase in chondrocytes, the inflammatory response, calcium metabolism and chondrocyte survival 32 . These observations suggest that MTCs are active in energy supply. In a previous study, signs of mitochondrial dysfunction were identified in cultured KBD chondrocytes, including decreased cellular ATP levels, an increased ratio of cells with de-energised mitochondria, release of mitochondrial cytochrome c and elevated activation of caspases 9 and 3. In addition, the number of apoptosis-positive chondrocytes in patients with KBD was larger than that in healthy controls 24,31 . All these findings suggest that mitochondrial dysfunction and mitochondriamediated cell death play critical roles in the pathophysiology of KBD.
We also identified a second novel cell population, HomCs. These cells were mainly present among KBD chondrocytes and exhibited high expression levels of metallothionein (MT) genes. MTs play diverse roles in metal detoxification, antioxidant functions, immune defence and cell differentiation. MTs are also important in regulating intracellular ROS levels 33,34 . A previous study also indicated that oxidative stress is the common pathological change in juvenile KBD patients recruited from both severely affected and mildly affected endemic regions. Severe oxidative stress is also associated with the occurrence and development of KBD 35 . The evidence above suggests that these HomCs might play a protective role against apoptosis by scavenging ROS and controlling cellular oxidative stress in KBD.
Although lacking intrinsic reparative ability, articular cartilage contains a population of cells with progenitor- like qualities [36][37][38][39][40][41] that are similar to stem cells. 42 . This population of chondrocytes is known as CPCs and has been observed and isolated in human articular cartilage and characterised by their capacity for self-renewal, expression of stem cell-related surface markers and ability to differentiate along multiple lineages 17 . Since Dowthwaite et al. 43 first reported the isolation of CPCs from the surface zone of articular cartilage, CPCs have been detected mostly in the superficial zone of articular cartilage in a number of studies 44,45 . CPCs have also been found to participate in cartilage self-repair during early or later stages of human OA 20,46,47 and were found to be highly migratory towards damaged cartilage tissue and repopulated in repair tissue. However, another study indicated that PRG4-expressing cells served as a progenitor population of cartilage and were mainly found in superficial chondrocytes in young mice but expanded into deeper regions of articular cartilage as animals aged 48 . Yu et al. 49 demonstrated that CPCs can be isolated from both the superficial and deep bovine articular cartilage and identified CPCs in healthy articular cartilage by a singlecell clonogenicity screening technique for the first time. In our study, the CPC population was mainly found in the middle and deep zones of normal healthy cartilage and possessed a high proportion of genes related to selfrenewal ability, differentiation multipotency and migration 20,50 . The expression of CPC markers (PTTG1 and BIRC5) in OA and KBD cartilage tissues showed that CPCs were distributed in the superficial zone in OA cartilage and in the deep zone in KBD cartilage. These inconsistent results for the expression of CPC markers in OA and KBD may be caused by their different injury patterns to cartilage; OA is characterised by progressive degradation on the cartilage surface, while the degenerative changes in KBD cartilage are characterised by multiple foci of chondronecroses in the deep zone of the cartilage 24 . On the basis of these data, we inferred that the CPCs from the superficial zone in OA might actually originate or migrate from the middle and deep zones of CPCs in healthy articular cartilage. Moreover, CPCs could migrate extensively towards damaged cartilage tissue parts and repopulate in repair tissue.
Because healthy controls and KBD patients now rarely undergo surgical cartilage biopsy, our study was limited by a relatively small sample size. Cartilage samples were collected from only three pairs of patients and healthy controls, which may lead to statistical bias caused by the individual characteristics of the limited number of subjects. In addition, there is no evidence of sex differences in KBD, but a large number of studies have shown that women have a higher prevalence of knee OA than men, particularly after the age of 50 51,52 . In this study, to avoid the statistical bias caused by sex differences, all three cartilage samples used in scRNA-seq analysis were collected from female donors; however, we cannot ignore the influence of the sex bias of knee OA on the outcomes of the scRNA-seq analysis in future analyses. Therefore, we performed IHC to verify the results of scRNA-seq analysis and confirm similar results in cartilage from both women and men.
In conclusion, our study systemically compared the single-cell landscapes of KBD, OA and healthy chondrocytes. This analysis identified discriminative markers for specific chondrocyte populations, which may be considered targets for cartilage injury repair or key molecules in the pathogeneses of these two diseases. A new population, MTCs, was markedly expanded in KBD, and this expansion may represent an important mechanism of chondrocyte degeneration in KBD.

Sample recruitment
The patients with KBD were diagnosed strictly according to the national diagnostic criteria of KBD in China [WS/T 207-2010]. The patients with OA were diagnosed strictly according to the Modified Outerbridge Classification. All subjects with alterations such as defects and sclerosis on the bone end of phalanges combined with compression changes of the calcaneus and talus and enlarged/deformed limb joints (hand, elbow, knee, ankle, etc.) manifested on X-ray were diagnosed with KBD. Subjects were excluded if they were suffering or had previously suffered from any other osteoarticular diseases (such as rheumatoid arthritis, gout or skeletal fluorosis) or any other type of macrosomia, disorder of osteochondrodysplasia or chronic disease (such as hypertension, diabetes or coronary heart disease) or had received any treatment in the past six months. Clinical information was collected from patient records. Articular cartilage samples from KBD and OA patients were collected from individuals who underwent arthroplasty of the knee. Healthy controls were obtained from the patients who suffered trauma or amputation due to an accident.
All donors signed a written informed consent form. All subjects were of Chinese Han lineage. Adult articular cartilage specimens were collected from one KBD patient, one OA patient and one healthy subject (Table S1) for single-cell RNA-seq and from five KBD patients, five OA patients and five healthy subjects (Table S2) for immunohistochemistry verification and from three KBD patients, three OA patients and three healthy subjects (Table S2) for qRT-PCR verification.

Cartilage tissue collection and chondrocyte isolation
All articular cartilage samples, including all of the cartilage zones (including calcified) and subchondral bone, were harvested from the lateral tibial plateau and obtained within 1 h after operation. Chondrocytes were isolated as follows: 53-56 articular cartilage specimens were washed twice with sterile phosphate-buffered saline (PBS) with antibiotics (penicillin and streptomycin), cut into pieces (1 mm 3 ), and subjected to enzymatic digestion with 0.25% trypsin at 37°C in an atmosphere of 5% CO 2 for up to 30 min. Cell suspensions were centrifuged at 1000 × g for 5 min, the supernatant was aspirated completely, and the cells were digested in basal media supplemented with 0.2% type II collagenase at 37°C using an Eppendorf Thermomixer for 12-16 h. Isolated chondrocytes were filtered through 70 mM nylon filters, washed twice with sterile PBS and then directly prepared for cDNA amplification and scRNA-Seq library construction.
Single-cell RNA sequencing Cell capture and cDNA synthesis The Chromium Single Cell 3′ Library and Gel Bead Kit V2 (10x Genomics, 120237) and Single Cell A Chip Kit (10x Genomics, 120236) were used for cell capture. The cell suspension (300-600 living cells per microlitre, as determined by Count Star) was loaded onto the Chromium Single Cell Controller (10x Genomics) to generate single-cell gel beads in emulsion (GEMs) according to the manufacturer's protocol. In short, single cells were suspended in PBS containing 0.04% BSA. Approximately 3173 cells were added to each channel, and the target cell recovery was estimated to be~16,544 cells in total. Captured cells were lysed, and the released RNA was barcoded through reverse transcription in individual GEMs. Reverse transcription was performed on a S1000TM Touch Thermal Cycler (Bio Rad) at 53°C for 45 min, followed by 85°C for 5 min and a hold at 4°C. cDNA was generated and then amplified, and quality was assessed using an Agilent 4200 (performed by CapitalBio Technology, Beijing).

Single-cell RNA-Seq library preparation
According to the manufacturer's instructions, singlecell RNA-seq libraries were constructed using the Single Cell 3′ Library Gel Bead Kit V2. The libraries were sequenced using an Illumina NovaSeq6000 sequencer with a sequencing depth of at least 57,374 reads per cell with a paired-end 150 bp (PE150) reading strategy (performed by CapitalBio Technology, Beijing).

Seurat pipeline
Clustering was also performed with Seurat 3.0 (R package). Cells with a gene number was <200, a gene number ranked in the top 1%, or a mitochondrial gene ratio of >25% were regarded as abnormal and filtered out. Dimensionality reduction was performed using principal component analysis (PCA), and visualisation was realised by TSNE and UMAP.

Enrichment analysis
GO enrichment, KEGG enrichment, Reactome enrichment and disease enrichment (human only) of cluster markers were performed using KOBAS software with Benjamini-Hochberg multiple testing adjustment, using the top 20 marker genes for each cluster. The results were visualised in R.

Single-cell trajectory analysis
Single-cell trajectories were built with Monocle (R package) for pseudotime analyses. Genes were filtered by the following criteria: expressed in more than 10 cells; average expression value >0.1; and Qval < 0.01 in different analyses.

WGCNA
Weighted correlation network analysis (WGCNA) was performed by the WGCNA R software package with default parameters. According to the above clustering results, every cluster was divided into subclusters, and the average expression of genes in a subcluster was calculated.

Cell type annotation
Cell type annotation was carried out by SingleR (https:// bioconductor.org/packages/devel/bioc/html/SingleR. html), which performs unbiased cell type recognition from single-cell RNA sequencing data by leveraging reference transcriptomic datasets of pure cell types to infer the cell of origin of each single cell independently. The human blueprint_encode and Human Primary Cell Atlas (HPCA) reference datasets were used.

RNA extraction and quantitative reverse-transcription PCR verification
Total RNA was isolated from chondrocytes using the Trizol protocol. The RNA was then converted into complementary DNA (cDNA) by RT-PCR using the RevertAid™ First Strand cDNA Synthesis Kit (Thermo Scientific Molecular Biology, Canada) according to the manufacturer's instructions. qRT-PCR was performed using an ABI7500 Real-Time PCR system (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's instructions. All primer and probe sets were supplied by TaqMan® Gene Expression Assays (Applied Biosystems). The relative expression level of the gene was calculated using the comparative Ct method. β-actin was used as an internal control to normalise the sample differences.

Immunohistochemical verification
Cartilage tissues were fixed with 4% (w/v) paraformaldehyde for 24 h immediately after acquisition and decalcified in 10% (w/v) ethylenediaminetetraacetic acid disodium salt (EDTA-Na 2 ) for 2-3 weeks. The samples were dehydrated in an alcohol series, cleared in xylene and embedded in paraffin wax. Paraffin sections were cut into 5-µm sections, mounted on slides and stored at room temperature until use for staining. The paraffinembedded sections were baked at 65°C for 1 h, deparaffinized with xylene and then rehydrated in decreasing concentrations of ethanol. Endogenous peroxidase activity was blocked by 0.3% (w/v) hydrogen peroxide for 10 min at room temperature, and then the sections were washed with 1×PBS. Then, the sections were incubated in 10 mol/ L urea solution diluted with water at 37°C for 20 min and washed with 1×PBS. For analysis of extracellular markers, such as IGFBP5 and CHI3L1, the sections were treated with 2 mg/ml hyaluronidase at pH 5.0 and 37°C for 30 min. For analysis of intracellular markers, the sections were incubated in 0.1% trypsin/CaCl 2 at 37°C and digested for another 20 min for antigen retrieval. antibodies and IgG as a negative control were applied on the sections, and the samples were further incubated overnight at 4°C. After washing with 1×PBS, sections were incubated using a human serum amyloid P (SAP) kit (solution B contains a goat anti-rabbit secondary antibody; Zhongshan, Jinqiao, Guangzhou, China) according to the manufacturer's instructions. The substrate 3,3′diaminobenzidine was added to stain the sections, and haematoxylin counterstaining was performed. Finally, the sections were dehydrated and mounted under alcoholcleaned coverslips. All IHC staining was assessed under light microscopy by two pathologists who were blinded to the origin of the samples. Each zone was identified based on previously reported characteristics 39,57-59 and articular cartilage was divided into three cell morphologies, namely, the superficial, middle and deep zones, according to light microscopy observation. Chondrocytes in the superficial zone were relatively small and flat and were oriented with the long axis parallel to the surface; cells in the middle zone were larger and more rounded and were randomly distributed in the matrix with fibres running in oblique directions; and cells in the deep zone were larger in size and were arranged in a columnar manner perpendicular to the surface. Assessment of the staining throughout each cartilage zone included systematic counting of positive and negative cells, starting from the cartilage surface and progressing down through all layers of cartilage. Five randomly chosen fields in each zone were counted at ×50 magnification. The percentage of positive cells was calculated using the number of positively stained cells divided by the total number of cells (positively and negatively stained cells) in the chosen fields of view. The percentages of positive cells in different zones were calculated for each case and then for the different groups.

Statistical analysis
Statistical analysis was performed using the SPSS 18.0 package. Differences in means were evaluated by one-way analysis of variance (ANOVA) with Tukey's post hoc test for multiple comparisons. Student's t test was applied to evaluate the difference between two groups. A normality test was applied first for continuous variables before any further comparison analyses. Nonparametric methods (e.g., Mann-Whitney U for pairwise comparisons) were used when the data were not normally distributed. A P value < 0.05 was considered to suggest a significant difference. The results of IHC and qRT-PCR are presented in bar charts using GraphPad PRISM 6 (GraphPad, San Diego, CA, USA).

Ethics statement
This study was approved by the ethics committee of Xi'an Jiaotong University, and the study was performed in accordance with the Declaration of Helsinki.

Informed consent
Informed written consent was obtained from each patients.