Single-cell RNA sequencing reveals different chondrocyte states in femoral cartilage between osteoarthritis and healthy individuals

Background: Cartilage injury is the main pathological manifestation of osteoarthritis (OA). Healthy chondrocyte is a prerequisite for cartilage regeneration and repair. Differences between healthy and OA chondrocyte types and the role these types play in cartilage regeneration and OA progression are unclear.

Method: This study conducted single-cell RNA sequencing (scRNA-seq) on the cartilage from normal distal femur of the knee (NC group) and OA femur (OA group) cartilage, the chondrocyte atlas was constructed, and the differences of cell subtypes between the two groups were compared. Pseudo-time and RNA velocity analysis were both performed to verify the possible differentiation sequence of cell subtypes. GO and KEGG pathway enrichment analysis were used to explore the potential functional characteristics of each cell subtype, and to predict the functional changes during cell differentiation. Differences in transcriptional regulation in subtypes were explored by single-cell regulatory network inference and clustering (SCENIC). The distribution of each cell subtype in cartilage tissue was identified by immunohistochemical staining (IHC).

Result: A total of 75,104 cells were included, they were divided into 19 clusters and annotated as 11 chondrocyte subtypes, including two new chondrocyte subtypes: METRNL+ and PRG4+ subtype. METRNL+ is in an early stage during chondrocyte differentiation, and RegC-B is in an intermediate state before chondrocyte dedifferentiation. With cell differentiation, cell subtypes shift from genetic expression to extracellular matrix adhesion and collagen remodeling, and signal pathways shift from HIF-1 to Hippo. The 11 subtypes were finally classified as intrinsic chondrocytes, effector chondrocytes, abnormally differentiated chondrocytes and dedifferentiated chondrocytes. IHC was used to verify the presence and distribution of each chondrocyte subtype.

Conclusion: This study screened two new chondrocyte subtypes, and a novel classification of each subtype was proposed. METRNL+ subtype is in an early stage during chondrocyte differentiation, and its transcriptomic characteristics and specific pathways provide a foundation for cartilage regeneration. EC-B, PRG4+ RegC-B, and FC are typical subtypes in the OA group, and the HippO-Taz pathway enriched by these cell subtypes may play a role in cartilage repair and OA progression. RegC-B is in the intermediate state before chondrocyte dedifferentiation, and its transcriptomic characteristics may provide a theoretical basis for intervening chondrocyte dedifferentiation.

1 Introduction

Osteoarthritis (OA) is one of the most common chronic diseases that cause disability in the elderly. More than 200 million people are affected by OA, which brings a great social burden (1). Traditionally, OA pathological changes are mainly the degeneration and loss of peri-articular cartilage and abnormal proliferation of subchondral bone (2). Therefore, clarifying the physiological and pathological characteristics of chondrocytes can provide the most direct reference for seeking ideal seed cells for cartilage regeneration (3). However, the role of OA chondrocytes in disease progression and the role of normal chondrocytes in cartilage regeneration have not been fully studied.

As a connective tissue, articular cartilage serves as lubricating and absorbing force via extracellular matrix (ECM), which accounts for about 95% volume of the cartilage, and is mainly composed of proteoglycans, collagen and water (4, 5). The cross-linking of proteoglycan and collagen is the key process to absorb force. Proteoglycan is thought to contribute more to absorbing force due to its binding to water (6), while collagen is considered to be more helpful in maintaining the integrity of cartilage tissue because of its cross-linking properties (7). Cartilage tissue can be divided into superficial, radial, deep, and calcified layers. In normal cartilage, the high percentage of collagen in the surface layer is beneficial to cope with the shear force, and the high content of proteoglycan in the radial layer and deep layer helps to cope with the pressure from the surface layer (8). The appearance of cartilage surface lesions in the early stage of OA indicates that collagen is one of the early stages of cartilage destruction (9). Account for around 5% volume of cartilage, chondrocytes maintain the turnover of ECM mainly through the balance of synthesis and catabolism. The imbalance of this process and the change of chondrocyte matrix gene expression profile are important factors causing cartilage tissue degeneration, and also the key targets in achieving cartilage regeneration (10–12). However, the functional heterogeneity of chondrocyte subtypes and its relationship to cellular synthesis and cartilage regeneration have not been adequately classified.

Single-cell RNA sequencing (scRNA-seq) characterizes the chondrocyte subtypes with different functional properties by clustering the transcriptomes of individual cells. Our previous scRNA-seq studies of cartilage from rheumatoid arthritic knee (13) and healthy ankle (14) identified two immune-associated chondrocyte subtypes and generated healthy ankle cartilage atlas, respectively. Tang et al. (15) conducted scRNA-seq on chondrocytes from OA cartilage and identified chondrocyte subsets, such as fibrochondrocytes (FC), effector chondrocytes (EC), homeostasis chondrocytes (HomC), regulatory chondrocytes (RegC), etc, for the first time. However, because of the small number of included chondrocytes, they did not analyze cartilage regeneration and OA related chondrocyte subtypes.

In order to screen chondrocyte subtypes related to cartilage regeneration and OA, we conducted scRNA-seq on a large number of chondrocytes from OA and normal femoral cartilage of the knee. This study aims to explore potential OA and cartilage regeneration associated chondrocyte subsets and supplement and verify previous studies. In addition, this study aims to deeply understand the differentiation sequence of each chondrocyte subtype, and explore the branching point of chondroblast subtype and abnormal differentiation subtype. Lastly, this study hopes to provide a theoretical basis for further exploring the pathogenesis of OA and seed cells for cartilage regeneration.

2 Methods

2.1 Volunteer screening and sample selection

Samples were obtained from OA patients who underwent TKA and normal knee joint population for corpse donation. The inclusion criteria for patients were as follows: a diagnosis of primary osteoarthritis of the knee with K-L grade III or IV. For knee joints donated from corpses, samples were evaluated as healthy according to the Modified Outerbridge Classification before being included in this study. Exclusion criteria included genu valgus or genu varus ≥ 20°, inflammatory arthritis; gouty arthritis, congenital knee dysplasia and a history of previous knee surgery. The full layer of cartilage deep to the bone surface was obtained for scRNA-seq. Finally, we obtained 3 NC samples and 3 OA samples respectively. Next, we performed cell isolation, cell counting and quality control, scRNA-seq and data analysis (Figure 1A). This study was approved by the institutional ethics review board at the Affiliated Hospital of Qingdao University and all donors signed written informed consent.

Figure 1

Figure 1 The process of scRNA-seq and the results of cell annotation (A) Flow chart of the experimental strategy. It includes sampling from femur of NC and OA, single cell isolation, cell counting and quality control, sequencing library preparation, single cell RNA-seq, phenotypic classification and statistical analysis. (B) Coefficient for nCount_RNA and nFeature_RNA (C) Parameters for screening nCount_RNA, nFeature_RNA and mitochondrial genes (D) Umap of all cells colored according to cell clusters. The annotation results of 19 clusters are listed on the right. (E) Umap of all cells colored according to cell subtypes (F) Heat map showing scaled differently expressed genes for each cell type defined in (E). Specific representative genes in each cell type are listed along the left of heat map. (G) umap of all cells according to indicated marker genes for each cell subtype.

2.2 Tissue dissociation and preparation of single-cell suspensions

Fresh cartilage tissue was transferred into a sterile culture dish containing phosphate-buffered saline (PBS) on ice, and cut it into 0.5 mm3 pieces. These pieces were washed with PBS twice and dissociated into single cells in dissociation solution (0.2% collagenase II 0.25% Trypsin-EDTA) in 37°C water bath with shaking for 20 min at 100 rpm. Digestion was terminated with PBS containing 10% fetal bovine serum (FBS).The resulting cell suspension was filtered by passing through 70–30 um stacked cell strainer and centrifuged at 300g for 5 min at 4°C. The cell pellet was added 100μl Dead Cell Removal MicroBeads (MACS 130–090-101) to remove dead cells using Miltenyi ® Dead Cell Removal Kit (MACS 130–090-101). Then, it was resuspended and certifuged twice at 300 g for 3 min at 4°C. The overall cell viability was confirmed by trypan blue exclusion, which needed to be above 85%. Single cell suspensions were counted using a Countess II Automated Cell Counter with a concentration of 700–1200 cells/μl.

2.3 10× Genomics scRNA-seq library preparation and sequencing

Single-cell suspensions were loaded to 10x Chromium to capture single cell according to the manufacturer’s instructions of 10x Genomics Chromium Single-Cell 3’ kit (V3). Steps of cDNA amplification and library construction were performed according to the standard protocol. Libraries were sequenced on an Illumina NovaSeq 6000 sequencing system (paired-end multiplexing run,150bp) by LC-Bio Technology co.ltd., (HangZhou, China) at a minimum depth of 20,000 reads per cell.

2.4 Bioinformatics analysis

Sequencing results were demultiplexed and converted to FASTQ format using Illumina bcl2fastq software (version2.20). We used Cell Ranger (https://support.10xgenomics.com/single-cell-geneexpression/software/pipelines/latest/what-is-cell-ranger, version 3.1.0) to analyze the quality of the original sequencing data of each sample and compared them with the reference genome (Ensembl genome GRCh38/GRCm38 reference genome). The analysis results of Cell Ranger were loaded into Seurat software (version3.1.1) for dimensional reduction and cluster analysis. Quality control was performed strictly. Firstly, the gene was detected in at least 3 cells were included. Secondly, the number of genes expressed in a single cell must be more than 500. Thirdly, the number of UMI must be greater than or equal to 500. Fourthly, the percent of mitochondrial-DNA derived gene-expression should be less than 25%.

To visualize the data, we further reduced the dimensionality of filtered cells using Seurat and used t-SNE to project the cells into 2D space. The steps included: using the Log Normalize method of the “Normalization” function of the Seurat package to calculated the expression value of genes; using the normalized expression value to perform PCA (Principal component analysis) analysis; using the top 10 PCs to do clustering and t-SNE analysis; selecting weighted Shared Nearest Neighbor (SNN) graph-based clustering method to find clusters. Marker genes for each cluster were identified with the Wilcoxon rank-sum test (default parameters is “bimod”: Likelihood-ratio test) via the FindMarkers function in Seurat. In addition, the gene expressed in more than 10% of the cells in a cluster and the gene expression multiple of log2FC ≥ 0.26 was identified as a marker gene. Differentially expressed genes in one cluster between HLB and LLB were also identified via the FindMarkers function in Seurat. Adjusted P value <0.05 and |log2 FC|>0.26 were set as cutoffs.

2.5 Gene ontology and Kyoto encyclopedia of genes and genomes pathway enrichment analysis

GO enrichment analysis provides significant GO terms by comparing differentially expressed genes with the whole gens of genome background (http://www.geneontology.org/). These significant GO terms (p-value < 0.05) provide clues for understanding the possible biological functions of cell population.

Genes usually interact with each other to play roles in certain biological functions. Pathway- based analysis helps to further understand biological functions of genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) is the major public pathway-related database. Pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways by comparing differentially expressed genes with the whole genome background genes.

2.6 Identification of different cell types

We further selected the top 10 genes as the marker genes according to differentially expression genes in each cluster. Then the distribution of each marker gene was demonstrated by using heat map, bubble diagram and umap. By comparing our marker genes and functional enrichments with that obtained from previous studies related to cartilage, we identified the cell type of each cluster.

2.7 Pseudotime analysis and RNA velocity analysis

Pseudotime analysis and trajectory construction were conducted by Monocle (http://cole-trapnell-lab.github.io/monocle-release/, version 2). The generated Seurat object was transformed into Monocle object with the “GetAssayData” and “newCellDataSet” function. The “orderCells” function in Monocle package arranged cells along a time axis to present their suggested position in a differential trajectory. The parameter of root_state was defined as ProC. The identified trajectory was colored by cell types or clusters to show the development pathway during chondrocyte differentiation.

2.8 Single-cell regulatory network inference and clustering

Seurat objects with annotations were imported to run Scenic (version 0.9.1). Use GENIE3 to extrapolate gene regulatory networks from gene expression data. The TF targeting relationship with the gene was verified with the help of RcisTarget software. Then, cells with active gene regulatory networks were identified using AUCell software, and regulons were scored for their Area Under Curve (AUC) by binarizing specific regulons. The RSS score of the TF was obtained by the calcRSS and ranked by value to determine the Top TF for each cell type. The differentially expressed transcription factors were visualized using the pheatmap function.

2.9 Immunohistochemical assays and cellular immunofluorescence detection

The immunohistochemical staining of trimmed articular cartilage tissue includes: fixing in 4% buffered paraformaldehyde for 48 hours, decalcifying with EDTA for about 4 weeks, embedding and sectioning in paraffin, dewaxing and rehydrating paraffin sections with xylene and gradient alcohol, blocking endogenous peroxidase activity with 3% H2O2 for 10 minutes, serum sealing with 3% BSA for 30 minutes at room temperature, antigen retrieval with pepsin digestion, primary antibody incubation overnight at 4°C, corresponding secondary antibody incubation for 50 minutes at room temperature, DAB color reaction and counterstaining nucleus with hematoxylin stain solution for about 3 minutes. Upright microscope (NIKON Eclipse ci) and imaging system (NIKON digital sight DS-FI2) were applied to enlarge and take photos for the histological section using white light. The widely accepted German semi-quantitative scoring system was used, taking into account the range of staining strength and staining area. Staining intensity: 0, no coloring; 1, light yellow; 2, brown-yellow; 3, brown- brown. Stained area: 0, < 5%; 1, 5% -25%; 2, 25%-50%; 3, 51%-75%; 4, >75%. These two scores were multiplied together as the final score. Comparative analysis of IHC staining scores for each layer of cartilage tissue was performed using the Kruskal-Wallis test for within-group differences and the Wilcoxon rank sum test for between-group differences.

The immortalized human chondrocytes (iCell Bioscience Inc, Shanghai) were passaged onto glass slides in a 6-well plate and grown to an experimental usable state. After washing with PBS buffer, the cells were fixed with 4% paraformaldehyde for 15 minutes. ANGPTL4 and CRLF1 primary antibodies were added separately and incubated overnight at 4°C to enhance antibody-cell interactions. After washing to remove unbound primary antibodies, corresponding secondary antibodies were added and incubated for 1 hour at 4°C. Following another wash to remove unbound secondary antibodies, anti-fade mounting medium was added to prevent signal fading. The slides were covered with glass coverslips, and fluorescence imaging was performed using a fluorescence microscope.

2.10 Statistical analysis

SPSS version 23.0 (IBM Inc., Armonk, NY, USA) was used to compare the difference in cell numbers between the NC group and OA group. T test was applied for testing, and P < 0.05 was considered to be statistically significant.

3 Results

3.1 Single-cell profiling and identification of femoral cartilage chondrocytes

75,104 cells were analyzed and grouped into clusters 0–18 (Figure 1D). The condition for quality control is shown in Figures 1B, C. The 19 clusters were then summarized and annotated into 11 distinct cell subtypes, including EC-A, EC-B, HTC-A, HTC-B, RegC-A, RegC-B, FC, METRNL+ subtype, PRG4+ subtype, C1QA+ subtype and an undefined subtype (Figure 1E). Figure 1F presents some highly expressed DEGs in each subtype. The marker genes for subtypes discussed in this study are shown in Figure 1G, which shows that CHRDL2 for EC-A, COL2A1 for EC-B, HILPDA for HTC-A, IBSP for HTC-B, VCAM1 for RegC-B, METRNL and ANGPTL4 for MERNL+subtype, PRG4 and CRLF1 for PRG4+ subtype.

Besides referring to marker genes reported in the literatures, we also considered the results of enrichment analysis (Figure 2A) to annotate cell subtypes. We summarized the top 20 terms ranked by p-value of each subtype in the biological process of GO enrichment analysis and presented them by heatmap. It was found that the subtypes that have a higher percentage in the NC group (Figures 2B, C) mainly presented hierarchical changes around GO terms related to the genetic expression process (Figure 2A). Although marker gene for the undefined subgroup was not screened successfully, the results of enrichment analysis showed that it was more prominent in the gene and chromatin damage repair. Although the functions including transcription, translation, protein modification, and transport for the remaining five subtypes were all significant, their specificities were different. METRNL+ subtype tended to play an important role in protein processing and transport, HTC-B was involved in endoplasmic reticulum translation, and HTC-A is preferentially enriched in transcription and ribosome processing. EC-A exhibited reactive activity, including biological processes such as response to hormones, regulation of the cell cycle, and differentiation into skeletal muscle cells. The GO enrichment results of RegC-A (cluster 11 and 14) are partly similar to the metal ion associated chondrocytes (MIRC) mentioned in the previous study of our team, which was enriched in response to metal ions (14). However, in combination with the differentiation sequence for RegC-A and RegC-B suggested by RNA velocity analysis shown in Figure 3F, and the specific expression of RegC marker genes such as CHI3L1, CHI3L2 and VCAM1 mentioned by Tang et al (15) in cluster11 and 14, we annotated them as RegC-A.

Figure 2

Figure 2 Enrichment analysis and proportion of each subtype, and RNA velocity analysis of two newly identified subtypes (A) GO enrichment analysis of each subtype. The top 20 biological process terms for each subgroup based on p-value were summarized and visualized by heatmap according to -log10p values. (B) The proportion of each subtype in each sample (C) t-test for the proportion of the number of subtypes in each sample. *p < 0.05, **p <0.01, ***p <0.001 (D) RNA velocity analysis of METRNL+ subtype, HTC-A, HTC-B, and EC-A. Arrows show the direction and rate of cellular movement. (E) RNA velocity analysis of EC-B, RegC-B, and PRG4+ subtypes. (F) The top 10 DEGs of METRNL+ subtype based on log2FC. (G) GO enrichment analysis of METRNL+ subtype. The gene with Log2FC*PCT1/PCT2≥1.5 in this subtype was applied to enrichment analysis. (H) The top 10 DEGs of PRG4+ subtype based on log2FC. (I) GO enrichment analysis of PRG4+ subtype. The gene with Log2FC*PCT1/PCT2≥1.5 in this subtype was applied to enrichment analysis.

Figure 3

Figure 3 Differentiation sequence of chondrocyte subtypes (A, D, G) pseudo-time analysis of all cells, NC group and OA group. The dots in each color represents each state, and each state is named numerically. The direction of the arrow represents the direction of the pseudo-time. (B, E, H) umap visualization of state in (A, D, G). The meanings of the colors, numbers, and arrows correspond to the (A, D, G), respectively. (C, F, I, J) RNA velocity analysis of EC-A, EC-B, and RegC-B in all cells, NC group, and OA group cells. (K) Relative expression levels of COL1/2 subtypes (L) Summary of the differentiation sequence of each subtype. The arrows represent the differentiation direction.

The subtypes that have a higher percentage in the OA group (Figures 2B, C) were mainly hierarchical enriched in GO terms related to damage repair (Figure 2A). EC-B was mainly enriched in responding to injury, such as ECM and collagen remodeling. Chondrocytes that specifically express COL1 have been considered as a sign of chondrocyte dedifferentiation, in the current study, FC owned this characteristic (Figure 3K). FC was preferentially enriched in cell migration and adhesion, ECM adhesion, and angiogenesis. RegC-B was mainly enriched in chondrocyte differentiation and development. Although the p-value was slightly lower, the GO terms of the PRG4+ subtype are similar to those of FC. In addition, we identified the GO items related to osteochondral biological function (Supplementary Figure 1A) with a p value of less than 0.05. The results showed that EC-B, RegC-B, PRG4+ and FC were more enriched in osteogenic function, cartilage differentiation and development than those in the subtypes that have a higher percentage in the NC group.

To compare the difference for each cell subtype between the NC and OA groups, we analyzed the proportion of cell subtypes in each sample (Figures 2B, C). The results showed that the proportion of EC-A, HTC-A, HTC-B, RegC-A and METRNL+ subtypes in the NC group was higher than that in the OA group (p<0.05), while the proportion of EC-B, RegC-B, FC and PRG4+ subtypes in NC group was lower than that in OA group (p<0.05). Considering the results of pseudo-time and RNA velocity analysis (Figure 3), we classified HTC-A, HTC-B, RegC-A, and METRNL+ subtype as inherent chondrocytes (InC), EC-A, and EC-B are classified as effector chondrocytes (EC), RegC-B and PRG4+ subtype as abnormally differentiated cells, and FC as dedifferentiated chondrocytes.

3.2 Biological characteristics METRNL+ and PRG4+ subtypes

As mentioned above, the number of METRNL+ subtype was much more in the NC group compared with the OA group, and its function such as the genetic expression process was particularly prominent. RNA velocity analysis of functionally similar subtypes (Figure 2D) showed that METRNL+ had the most active RNA metabolism, suggesting that it may be in the early stage during the differentiation process. Further exploration of the genetic expressional characteristics of METRNL+ subtype showed that besides METRNL, genes such as ANGPTL4, MMP3, LOX, PYGL, VEGFA, etc were also specifically highly expressed (Figure 2F, Supplementary Figure 2A). To explore the functional characteristics of METRNL+ subtype, DEGs of this subtype were further screened for more specific genes based on the results of Log2FC*PCT1/PCT2≥1.5 and then transferred to GO enrichment analysis (Figure 2G). The results suggested that this subtype is preferentially enriched in responding to hypoxia, angiogenesis, protein assembly, and differentiation into brown fat cells. In addition, signal transduction and apoptosis may be active in this subtype. We sorted out GO terms associated with bone and cartilage biological functions (Supplementary Figure 1A) and GO terms related to the METRNL gene (Supplementary Figure 1B) with a subset having a p-value less than 0.05. Our results demonstrate that, although this subtype exhibits enrichment in osteogenesis in METRNL+, the EC-B, RegC-B, PRG4+, and FC subgroups show a more pronounced enrichment in osteogenic functional terms based on p-value and rich factor. Notably, the METRNL+ subgroup displays a significant enrichment in chondroitin sulfate synthesis function. In summary, this subtype may play an important role in cartilage tissue maintenance and renewal.

We applied the same method to the PRG4+ subtype which has a lower proportion in the NC group compared with the OA group. RNA velocity analysis (Figure 2E) showed that RNA metabolism of PRG4+ subtype was slower than that of EC-B and RegC-B, suggesting that PRG4+ should be in a late stage during differentiation. PRG4+ subtype was preferentially enriched in GO terms related to cellular and tissue adhesion (Figure 2I), 24 genes such as FN1 and TNFAIP6 with Log2FC*PCT1/PCT2≥1.5 (Figure 2H, Supplementary Figure 2C) were associated with these biological processes. In addition, 15 genes, including PRG4, CRIP1 and CD55, were enriched in GO terms related to inflammation (Figure 2I). Furthermore, this subtype was also associated with ECM ossification, fibrogenesis, and TGFβ process (Figure 2I). In summary, this subtype may respond to injury and inflammation, and own ECM fibrogenesis and osteogenesis potential.

3.3 Pseudo-time and RNA velocity reveal the differentiation order of cell subtypes

We performed pseudo-time (Figures 3A, D, J) and RNA velocity analysis (Figures 3C, F, I) on all cells, NC group and OA group, respectively, and the state determined by pseudo-time analysis was further visualized by umap plot (Figures 3B, E, H). The differentiation sequence of each chondrocyte subtype is presented in Figure 2K. In detail, for all cells, cells in state1(METRNL+ chondrocyte, HTC, and EC-A) differentiated to both fate1 and fate2 directions, which represented state2 (undefined subtype) and state5 (EC-B, RegC-B, FC, PRG4+ chondrocyte), respectively (Figures 3A–C). This suggested a trend that the cells in the NC group differentiated into the OA group. The results of RNA velocity analysis of all cells also corresponded to this trend (Figure 3C).

Originating from state 6 (METRNL+ chondrocyte), cells in the NC group also differentiated to fate 3 (state 2, EC, and HTC) besides fate 1(Figures 3D–F). Consistent with the results of pseudo-time analysis, the RNA velocity analysis is presented in Figures 2D and 3F showed that the METRNL+ subtype was in an early stage during differentiation and differentiated to EC-A, RegC-B, and undefined subtype, respectively. As shown in Figures 3G–I, cells in the OA group originated from state8 (HTC-B and EC-B) and were further differentiated to FC via RegC-B. In combination with the results shown in Figure 2E, RegC-B also showed a tendency to differentiate to the PRG4+ subtype. In order to further clarify the differentiation relationship among EC-A, EC-B, and RegC-B, we conducted RNA velocity analysis on these cells (Figure 3J), and the results suggested that EC-A and EC-B could both differentiate to RegC-B, and EC-A tended to differentiate to EC-B.

In summary (Figure 3L), first, the METRNL+ subtype was in an early stage during differentiation, and the METRNL+ subtype in the NC group mainly differentiated into three directions: EC-A, RegC-B and undefined subtype. Second, FC was in a late stage that exhibited a dedifferentiated character, RegC-B was the major origin of FC, and EC-B was the main origin of RegC-B. In addition, EC-A and HTC-B were likely to be the main origins of EC-B. Last, part of EC-B and RegC-B were likely to differentiate to the PRG4+ subtype. We used the relative expression levels of COL1/COL2 to determine the dedifferentiation degree of each subtype (Figure 3K), which corresponded to the results of pseudo-time and RNA velocity analysis.

3.4 Functional changes during cell differentiation

To explore the functional evolution from intrinsic chondrocytes to abnormally differentiated chondrocytes, we performed transcriptomic comparisons and enrichment analysis on METRNL+ subtype, EC-A, EC-B, and RegC-B (Figure 4). We summarized the top 20 enrichment terms based on the p-value of biological processes in the GO enrichment analysis (Figures 4A, B), which showed that the METRNL+ subtype focused on RNA splicing, protein translation and transport, and apoptosis. EC-A preferred to participate in peptide hormone response, skeletal muscle cell differentiation, circadian rhythm, cytoplasmic translation initiation, and transcriptional regulation. EC-B were enriched in ECM (including collagen) remodeling and ossification, chondrocyte coagulation, and cell differentiation. RegC-B were preferentially enriched in the TGF-β pathway, blood vessel formation and development, wound healing, and cell adhesion.

Figure 4

Figure 4 Changes in enrichment analysis during cell differentiation (A, B) GO enrichment analysis of each subtype. The top 20 biological process terms based on p-value for each subtype were summarized and visualized by heatmap according to either the -log10p value (A) or the -log10Rich factor (B). (C–F) KEGG pathway enrichment analysis of METRNL+, EC-A, EC-B and RegC-B. The bottom horizontal axis represents the number of genes enriched into the each pathway. (G) Expression of important downstream target genes of Hippo signaling pathway in each subtype. The size of the circle represents the proportion of cells in the cell subtype expressing the gene, and the color of the circle represents the average gene expression. (H) Expression of WWTR1 in each sample.

The results of KEGG analysis suggested that the METRNL+ subtype was enriched in HIF-1 and Ras signaling pathway, ribosome biosynthesis, amino sugar and nucleotide sugar metabolism, and glycolysis (Figure 4C). EC-A (Figure 4D) were involved in cAMP and Apelin signaling pathways, RNA degradation and degradation of valine, leucine, and isoleucine, and oxidative phosphorylation. EC-B preferred to participate in cytoskeletal regulation, Hippo, TGF-beta, and ECM-receptor interaction signaling pathways, RNA polymerase and proteasome, oxidative phosphorylation, and purine metabolism-related items (Figure 4E). To explore the targets of the Hippo signaling pathway, we compared the expression of all downstream target genes among all subtypes (Figure 4G) and all samples (Figure 4H), respectively, and found that the TAZ target (WWTR1) was more highly expressed in EC-B, RegC-B, PRG4+, FC, and OA samples. Therefore, we speculated that HiPO-TAZ should be one of the important ways that induce or against OA. Unlike EC-B, RegC-B was also enriched in the Rap1 signaling pathway, pyrimidine and triglyceride metabolism, drug metabolism and other related terms (Figure 4F).

In general, with the process of METRNL+ differentiation to RegC-B, the priority of protein synthesis and apoptosis gradually decreased. However, the functions of collagen synthesis and ECM ossification, chondrocyte coagulation and differentiation, cell adhesion and wound healing, angiogenesis and development have gradually increased in priority.

3.5 Differences in transcriptional regulation during cell differentiation

To clarify the evolution of transcriptional regulation from EC to abnormally differentiated chondrocytes, we performed transcriptional comparisons and functional enrichment analyses on EC-A, EC-B, and RegC-B (Figure 5). Figure 5A shows the top 10 specific transcription factors for each chondrocyte subtype. After RSS scoring, we obtained 5 prominent transcription factors for each subtype. Three significant transcription factors, MSX2 (EC-B), TCF4 (RegC-B), and DDIT3 (EC-A), were finally screened according to the expression levels of transcription factors based on umap plot (Figures 5C–E). To clarify the regulatory characteristics of the three transcription factors, we mapped the regulatory networks among the three transcription factors and DEGs of their subtypes.

Figure 5

Figure 5 Transcriptional regulation of chondrocyte subtypes (A) Heatmap of specific transcriptional factors of EC-A, EC-B, and RegC-B (B) RSS scores for specific transcription factors of EC-A, EC-B and RegC-B. (C–E) umap plot for key transcription factors of EC-A, EC-B, and RegC-B. (F) Gene regulatory network of specific transcription factors. The purple genes are the screened key transcription factor. The green and yellow genes are the highly expressed DEGs in the chondrocyte subtypes, and the yellow genes means that the gene itself is a transcription factor.

3.6 The distribution of subtypes in cartilage tissue

Similar to the percentage of chondrocyte subtypes in the two groups of scRNA-seq, IHC also indicated that the number of EC-A, HTC-A, and METRNL+ chondrocytes in the NC group was more than that in the OA group (p<0.05), and the number of RegC-B and FC in OA group was higher than that in OA group (p<0.05) (Figure 6). In addition, compared with the NC group, the number of EC-B and PRG4+ chondrocytes were more only in the middle layer of cartilage from the OA group (p<0.05) (Figure 6). Moreover, compared with the intermediate and deep layers, EC-A, EC-B, FC, METRML+, and PRG4+ chondrocytes were more abundant in the shallow layer (Figure 6). To further validate the two newly proposed subtypes, we conducted immunofluorescence staining on human chondrocytes. The results revealed widespread expression of the differentially expressed genes in these subtypes within human chondrocytes (Supplementary Figures 2B, C).

Figure 6

Figure 6 IHC results for specific markers for each cell type in different groups and layers. Scale bar, left, 200 μm; right, 50 μm. The scores of the indicated genes in the cartilage of both regions based on the immunohistochemistry assay are shown. *p<0.05, **p<0.01, ***p<0.001; ****p<0.0001; otherwise, not significant NC, normal femur cartilage; OA, OA femur cartilage.

4 Discussion

Chondrocytes are traditionally considered to have a poor regenerative ability and cannot be satisfactorily repaired after lesion. Recent studies tended to use stem cells and autologous chondrocytes in regenerating cartilage, but their efficacy remains to be further verified (16). This study screened METRNL+ cell subsets in the femoral cartilage of the knee for the first time. Although the role of METRNL in osteoblast function has been previously reported (17, 18), our results suggest that the osteogenic function is not the predominant function in the METRNL+ subtype. Its GO terms focus on functions such as protein synthesis and transportation, adding that representative genes of catabolism such as MMP3 are specifically expressed in this subtype, indicating that this subtype is a group of chondrocytes with active synthesis and catabolism process. The balance between synthesis and catabolism in chondrocytes is the driving force for cartilage tissue to maintain and renew its ECM (19, 20). Moreover, previous research has demonstrated the potential of METRNL in mitigating inflammation and pyroptosis in osteoarthritis (21), combined with METRNL+ subtype mainly comes from healthy cartilage, which suggest the cartilage protective effect of this subtype. Furthermore, pseudo-time analysis and RNA velocity analysis suggested that METRNL+ subtype belonged to early differentiated chondrocytes. In summary, the discovery of this subtype and its transcriptomic properties is expected to provide a new theoretical basis and practical direction for cell transplantation and cartilage regeneration.

The study conducted by Tang et al (15). mentioned that RegC might be involved in immune inflammatory response. A current study found that besides RegC, a group of chondrocytes with high expression of CRLF1, CD55, CRIP1, and PRG4 were also involved in the inflammatory response, and GO enrichment analysis showed that they tended to respond to acute inflammation. GO terms such as adhesion function of cells and extracellular matrix, osteogenesis and ossification regulation were also clearly significant in this subtype, which are closely related to the abnormal repair of cartilage injury (22–24). Therefore, it is expected to provide a new direction for preventing diseases including cartilage injury and OA by intervening in the abnormal repair of this subtype and verifying the mechanism for abnormal repair.

It is believed that EC is in the early stage during chondrogenic differentiation, and RegC is distributed throughout the chondrogenic differentiation process (15, 25). In this study, through pseudo-time and RNA velocity analysis, we suggested that the METRNL+ chondrocyte subtype should be distributed in a much earlier stage. EC owns a property of active compensatory, and its differentiation sequence may be after METRNL+ chondrocytes. RegC-B, as an intermediate state in the transition from nearly all OA and NC chondrocytes to FC, may be in an abnormal differentiation stage during the dedifferentiation process. COL1/2 is considered to be an important indicator of chondrocyte dedifferentiation (26), and the dedifferentiating degree of RegC-B in this study was between FC and EC. Therefore, the transcriptional regulation mechanism of RegC-B may provide a new direction for chondrocyte dedifferentiation.

Submerged in the joint fluid, chondrocytes are in a state of gradient hypoxia, and it has been found that the activation of the HIF-1 pathway, which is associated with hypoxia, stimulates the anabolism of chondrocytes (27, 28). In this study, the METRNL+ subtype was enriched in responding to hypoxia according to GO enrichment, and the HIF-1 pathway was also enriched in KEGG pathway, which further suggested the protective role of METRNL+ chondrocytes in the maintenance and regeneration of cartilage tissue. In addition, with the process of abnormally differentiation and dedifferentiation, the significant enriched pathway for RegC-B was transformed into the Hippo signaling pathway. Most scholars believe that the Hippo-YAP signaling pathway plays a protective role against (29, 30), but lack of reports on the role of the Hippo-TAZ pathway in the process of OA cartilage (31). The results of the current study suggest that Hippo-TAZ should be an important pathway to cause or against OA. Therefore, the mechanism Hippo-TAZ pathway on cartilage should be further explored.

Compared with previous studies, the current study conducted scRNA-seq on more chondrocytes and focused on the distal femur cartilage. In addition, the current study deeply analyzed each cell subtype in a customized way, speculated and classified the intrinsic chondrocytes, effector chondrocytes, abnormally differentiated chondrocytes, and dedifferentiated chondrocytes. Moreover, the current study explored the possible differentiation sequence and functional changing of the chondrocytes through pseudo-time and RNA velocity analysis, and proposed the potential state evolution of chondrocytes. Lastly, the typical cells of the NC and the OA groups and their transcriptomic characteristics were obtained via comparing the differences between the two groups, which is expected to provide potential targets for cartilage disease intervention and cartilage regeneration. The obvious differences between the two groups also provided straightforward results for future studies. However, the sample of mild cartilage lesion from OA was not included in this study, and the cellular atlas and transcriptomic characteristics of early OA may not be fully summarized.

5 Conclusion

METRNL+ is a typical subtype in the NC group, which is an early differentiation chondrocyte and may have a protective effect on cartilage tissue renewal and regeneration. EC-B, FC, PRG4+ and RegC-B cell subsets are typical cells in the OA group. The cartilage repair mediated by these subtypes is different from the renewal mechanism of ECM by normal chondrocytes, which may be the intrinsic factor of OA cartilage tissue degeneration. RegC-B is in the intermediate state before chondrocyte dedifferentiation, and its transcriptomic characteristics provide a theoretical basis for intervening chondrocyte dedifferentiation. HIF-1 is a typical pathway in normal chondrocytes, and the HiPO-TAz pathway may play a role in OA cartilage repair and OA progression.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: PRJNA1105421(SRA).

Ethics statement

The studies involving humans were approved by the Ethics Review Committee of the Affiliated Hospital of Qingdao University. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author contributions

ZS: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Software, Supervision, Validation, Writing – original draft, Writing – review & editing. MY: Conceptualization, Data curation, Investigation, Software, Methodology, Writing – review & editing. JW: Data curation, Methodology, Software, Conceptualization, Investigation, Writing – review & editing. HZ: Data curation, Methodology, Software, Writing – review & editing. XJ: Data curation, Methodology, Writing – review & editing. YX: Data curation, Methodology, Writing – review & editing. TW: Conceptualization, Investigation, Project administration, Resources, Supervision, Validation, Visualization, Writing – review & editing. TY: Conceptualization, Investigation, Project administration, Resources, Supervision, Validation, Visualization, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Acknowledgments

Thanking the participants for providing cartilage samples.

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/fimmu.2024.1407679/full#supplementary-material

Supplementary Figure 1 | GO terms related to osteochondral biological functions and the METRNL gene. (A) GO terms related to osteochondral biological functions. The color gradient represents the range of p-values from 0 to 0.01, while GO terms with p-values ranging from 0.01 to 0.05 are uniformly indicated in gray (B) GO terms related to the METRNL in METRNL+ subtype.

Supplementary Figure 2 | Immunofluorescence staining on human chondrocytes. (A, C) The top 10 DEGs of METRNL+ and PRG4+ subtype respectively based on log2FC. (B, D) Immunofluorescence staining of METRNL+ and PRG4+ subtype respectively. The three panels in the figure represent the immunofluorescence of antibodies (left panel), DAPI staining (middle panel), and the merged image of both (right panel), respectively. The images were captured at a magnification of 40X.

References

1. Steinmetz JD, Culbreth GT, Haile LM, Rafferty Q, Lo J, Fukutaki KG, et al. Global, regional, and national burden of osteoarthritis, 1990–2020 and projections to 2050: a systematic analysis for the Global Burden of Disease Study 2021, The Lancet. Rheumatology. (2023) 5:e508–22. doi: 10.1016/S2665-9913(23)00163-7

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Varela-Eirin M, Loureiro J, Fonseca E, Corrochano S, Caeiro JR, Collado M, et al. Cartilage regeneration and ageing: Targeting cellular plasticity in osteoarthritis. Ageing Res Rev. (2018) 42:56–71. doi: 10.1016/j.arr.2017.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Giorgino R, Albano D, Fusco S, Peretti GM, Mangiavini L, Messina C. Knee osteoarthritis: epidemiology, pathogenesis, and mesenchymal stem cells: what else is new? An Update. Int J Mol Sci. (2023) 24(7):6405. doi: 10.3390/ijms24076405

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Goldring SR, Goldring MB. Changes in the osteochondral unit during osteoarthritis: structure, function and cartilage-bone crosstalk, Nature reviews. Rheumatology. (2016) 12:632–44. doi: 10.1038/nrrheum.2016.148

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Krishnan Y, Grodzinsky AJ. Cartilage diseases. Matrix Biology: J Int Soc Matrix Biol. (2018) 71–72:51–69. doi: 10.1016/j.matbio.2018.05.005

CrossRef Full Text | Google Scholar

6. Nia HT, Han L, Bozchalooi IS, Roughley P, Youcef-Toumi K, Grodzinsky AJ, et al. Aggrecan nanoscale solid-fluid interactions are a primary determinant of cartilage dynamic mechanical properties. ACS Nano. (2015) 9:2614–25. doi: 10.1021/nn5062707

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Bielajew BJ, Hu JC, Athanasiou KA. Collagen: quantification, biomechanics, and role of minor subtypes in cartilage, Nature reviews. Materials. (2020) 5:730–47. doi: 10.1038/s41578-020-0213-1

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Ansari S, Khorshidi S, Karkhaneh A. Engineering of gradient osteochondral tissue: From nature to lab. Acta Biomaterialia. (2019) 87:41–54. doi: 10.1016/j.actbio.2019.01.071

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hosseini SM, Veldink MB, Ito K, van Donkelaar CC. Is collagen fiber damage the cause of early softening in articular cartilage? Osteoarthritis Cartilage. (2013) 21:136–43. doi: 10.1016/j.joca.2012.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Pap T, Bertrand J. Syndecans in cartilage breakdown and synovial inflammation, Nature reviews. Rheumatology. (2013) 9:43–55. doi: 10.1038/nrrheum.2012.178

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Cho Y, Jeong S, Kim H, Kang D, Lee J, Kang SB, et al. Disease-modifying therapeutic strategies in osteoarthritis: current status and future directions. Exp Mol Med. (2021) 53:1689–96. doi: 10.1038/s12276-021-00710-y

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Guilak F, Nims RJ, Dicks A, Wu CL, Meulenbelt I. Osteoarthritis as a disease of the cartilage pericellular matrix. Matrix Biol.. (2018) 71–72:40–50. doi: 10.1016/j.matbio.2018.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yan M, Sun Z, Wang J, Zhao H, Yu T, Zhang Y, et al. Single-cell RNA sequencing reveals distinct chondrocyte states in femoral cartilage under weight-bearing load in Rheumatoid arthritis. Front Immunol. (2023) 14:1247355. doi: 10.3389/fimmu.2023.1247355

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Wang T, Wang J, Sun Z, Zhang L, Yu C, Zhao H. Single-cell RNA sequence presents atlas analysis for chondrocytes in the talus and reveals the potential mechanism in coping with mechanical stress. Front Cell Dev Biol. (2022) 10:1047119. doi: 10.3389/fcell.2022.1047119

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ji Q, Zheng Y, Zhang G, Hu Y, Fan X, Hou Y, et al. Single-cell RNA-seq analysis reveals the progression of human osteoarthritis. Ann Rheumatic Dis. (2019) 78:100–10. doi: 10.1136/annrheumdis-2017-212863

CrossRef Full Text | Google Scholar

16. Campbell TM, Trudel G. Protecting the regenerative environment: selecting the optimal delivery vehicle for cartilage repair-a narrative review. Front Bioengineering Biotechnol. (2024) 12:1283752. doi: 10.3389/fbioe.2024.1283752

CrossRef Full Text | Google Scholar

17. Huang R, Balu AR, Molitoris KH, White JP, Robling AG, Ayturk UM, et al. The role of Meteorin-like in skeletal development and bone fracture healing. J orthopaedic research: Off Publ Orthopaedic Res Soc. (2022) 40:2510–21. doi: 10.1002/jor.25286

CrossRef Full Text | Google Scholar

18. Baht GS, Bareja A, Lee DE, Rao RR, Huang R, Huebner JL, et al. Meteorin-like facilitates skeletal muscle repair through a Stat3/IGF-1 mechanism. Nat Metab. (2020) 2:278–89. doi: 10.1038/s42255-020-0184-y

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Zelinka A, Roelofs AJ, Kandel RA, De Bari C. Cellular therapy and tissue engineering for cartilage repair. Osteoarthritis Cartilage. (2022) 30:1547–60. doi: 10.1016/j.joca.2022.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Chen Y, Cheng RJ, Wu Y, Huang D, Li Y, Liu Y. Advances in Stem Cell-Based Therapies in the Treatment of Osteoarthritis. Int J Mol Sci. (2023) 25(1):394. doi: 10.3390/ijms25010394

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Liu J, Jia S, Yang Y, Piao L, Wang Z, Jin Z, et al. Exercise induced meteorin-like protects chondrocytes against inflammation and pyroptosis in osteoarthritis by inhibiting PI3K/Akt/NF-κB and NLRP3/caspase-1/GSDMD signaling. Biomedicine Pharmacotherapy = Biomedecine Pharmacotherapie. (2023) 158:114118. doi: 10.1016/j.biopha.2022.114118

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Vonk LA, de Windt TS, Kragten AH, Beekhuizen M, Mastbergen SC, Dhert WJ, et al. Enhanced cell-induced articular cartilage regeneration by chondrons; the influence of joint damage and harvest site. Osteoarthritis Cartilage. (2014) 22:1910–7. doi: 10.1016/j.joca.2014.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Chen S, He T, Zhong Y, Chen M, Yao Q, Chen D, et al. Roles of focal adhesion proteins in skeleton and diseases. Acta Pharm Sinica. B. (2023) 13:998–1013. doi: 10.1016/j.apsb.2022.09.020

CrossRef Full Text | Google Scholar

24. Chen Y, Mehmood K, Chang YF, Tang Z, Li Y, Zhang H. The molecular mechanisms of glycosaminoglycan biosynthesis regulating chondrogenesis and endochondral ossification. Life Sci. (2023) 335:122243. doi: 10.1016/j.lfs.2023.122243

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Gao H, Di J, Yin M, He T, Wu D, Chen Z, et al. Identification of chondrocyte subpopulations in osteoarthritis using single-cell sequencing analysis. Gene. (2023) 852:147063. doi: 10.1016/j.gene.2022.147063

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Charlier E, Deroyer C, Ciregia F, Malaise O, Neuville S, Plener Z, et al. Chondrocyte dedifferentiation and osteoarthritis (OA). Biochem Pharmacol. (2019) 165:49–65. doi: 10.1016/j.bcp.2019.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Zhang H, Wang L, Cui J, Wang S, Han Y, Shao H, et al. Maintaining hypoxia environment of subchondral bone alleviates osteoarthritis progression. Sci Adv. (2023) 9:eabo7868. doi: 10.1126/sciadv.abo7868

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Stegen S, Laperre K, Eelen G, Rinaldi G, Fraisl P, Torrekens S, et al. HIF-1α metabolically controls collagen synthesis and modification in chondrocytes. Nature. (2019) 565:511–5. doi: 10.1038/s41586-019-0874-3

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Ma T, Ruan H, Lv L, Wei C, Yu Y, Jia L, et al. Oleanolic acid, a small-molecule natural product, inhibits ECM degeneration in osteoarthritis by regulating the Hippo/YAP and Wnt/β-catenin pathways. Food Funct. (2023) 14:9999–10013. doi: 10.1039/D3FO01902K

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Zhang S, Zhang B, Liao Z, Chen Y, Guo W, Wu J, et al. Hnrnpk Protects Against Osteoarthritis through Targeting WWC1 mRNA and Inhibiting Hippo Signaling Pathway. Mol Ther. (2024) 32(5):1461–78. doi: 10.1016/j.ymthe.2024.02.027

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Sun K, Guo J, Guo Z, Hou L, Liu H, Hou Y, et al. The roles of the Hippo-YAP signalling pathway in Cartilage and Osteoarthritis. Ageing Res Rev. (2023) 90:102015. doi: 10.1016/j.arr.2023.102015

PubMed Abstract | CrossRef Full Text | Google Scholar