Single-Cell RNA Sequencing Reveals Transcriptional Changes in the Cartilage of Subchondral Insufficiency Fracture of the Knee

Introduction

Spontaneous osteonecrosis of the knee (SONK) is a common disease that causes a sudden onset of knee pain. It was first described in the literature by Ahlback in 1968.¹ Later, SONK has been associated with SIFK once the subchondral fragment undergoes true necrosis; therefore, some authors believe the term “SONK” is a misnomer. Until now, the etiology of SIFK has been largely unknown. Low bone mineral density, varus deformity at the proximal tibia, and meniscal dysfunction were reported to be highly associated with the occurrence of SIFK.² For instance, medial meniscus posterior root tears (MMPRTs) and/or medial meniscus extrusion have been found to be highly correlated with SIFK.³ Recent studies have postulated that subchondral insufficiency fractures resulting from mechanical overloading concentrated on the medial part of the knee could cause SIFK.² This is supported by the histological finding that a fracture line occurs within the subchondral bone.⁴ Therefore, a more appropriate term for the condition may be subchondral insufficiency fracture of the knee (SIFK), rather than SONK.

Abnormalities of subchondral bone have been demonstrated to play a fundamental role in mediating cartilage degeneration in osteoarthritis (OA), which can be secondary to SIFK.⁵ In addition, changes in the medial meniscus (tear and/or extrusion) and/or subchondral bone insufficiency fracture might result in massive stress loading on the cartilage. Therefore, a previous study emphasized the need for assessment of the cartilaginous coating by means of arthroscopy before instituting treatment for SIFK.⁶ In fact, several studies have reported that the proliferative and redifferentiation abilities of chondrocytes from patients with SIFK are largely compromised compared with those of normal cells.^7,8 In agreement with this, levels of chondroitin 6-sulfate (C6S), which originates mainly from the articular cartilage, are found to be markedly elevated in the synovial fluid of SIFK patients compared to patients with medial OA, indicating cartilage breakdown.⁹ However, how a sudden insufficiency fracture of subchondral bone may contribute to cartilage metabolism alteration and/or ultimately lead to secondary OA is largely unclear. In this work, we were able to identify a subset of cell populations that were highly enriched with hypoxia and endochondral ossification molecules in SIFK cartilage compared with normal and OA cartilage. Furthermore, chondrocyte trajectory analysis and immunohistochemistry (IHC) staining demonstrated that transcription factors, including TCF4, were elevated in SIFK.

Materials and Methods

Clinical Specimen Collection

The diagnosis of SIFK was confirmed on MRI images when the T1-weighted image showed a low signal in the weight-bearing area and the T2-weighted image showed a high signal area surrounding the same site. The cartilage specimens were harvested from a 62-year-old female who underwent unicompartmental knee arthroplasty (UKA) surgery due to SIFK (Koshino classification: stage 2) and they were immediately subjected to chondrocyte isolation for scRNA-seq.

Chondrocyte Isolation and scRNA-Seq

The whole layers of SIFK cartilage specimens were subjected to collagenase type I (LS004196, 2mg/mL, 4h, 37°C) digestion to isolate the chondrocytes. Then, approximately 20,000 cells were used to acquire a single library by a Chromium controller (V3 chemistry version, 10X Genomics Inc, San Francisco, USA); the cells were barcoded and purified as we previously described,¹⁰ and then sequenced using a 2×150 paired-end configurations on the Illumina HiSeq platform at a sequencing depth of ~ 400 million reads. Cell Ranger (Version 6.0.2) was used to demultiplex the reads, which was followed by extraction of the cell barcodes and unique molecular identifiers (UMIs). The cDNA insert sequence was aligned to the human reference genome (GRCh38-2020).

Raw scRNA-seq files of normal and OA cartilage tissues were acquired from GEO: GSE169454.¹¹ Seurat package V3¹² was used for individual or integrated analysis of the datasets. Standard Seurat pipeline was used for filtering, variable gene selection, dimensionality reduction analysis and clustering. Doublets or cells with poor quality (genes >6000, genes <500, or >25% genes mapping to mitochondrial genome) were excluded. Expression was natural log transformed and normalized for scaling the sequencing depth to a total of 1×10⁴ molecules per cell. For the integrated dataset, batch integration was performed using Harmony (version 1.0).¹³ Briefly, the Seurat objects were first randomly downsampled to get comparable cell number between different group, then the objects were passed on to harmony integration by “RunHarmony” function, followed by downstream dimensional reduction by using the corrected Harmony embeddings. Only chondrocyte subpopulation remained for further analysis. Label transfer was done with Seurat “TransferData” function to classify the query SIFK chondrocyte cells based on reference OA chondrocyte (GSE104782) data,¹⁴ since the chondrocyte subtype of GSE169454 was also based on the clustering information of GSE104782. Gene set variation analysis (GSVA) was performed using the GSVA package for analyzing Hallmark and KEGG pathway signaling.¹⁵ Gene ontology (GO) analysis was performed using the clusterProfiler package.¹⁶ The cell metabolism state was inferred with the scMetabolism package,¹⁷ and a list of metabolic gene sets was acquired from the same package. Differentially expressed genes between normal and radiated groups were identified using FindMarkers within Seurat and p values between them were calculated by “test.use=bimod”.

Trajectory analysis was performed using Slingshot¹⁸ to computationally delineate the differential progression from normal cartilage to SIFK or OA cartilage. Briefly, Seurat objects were transformed into SingleCellExperiment objects. Slingshot trajectory analysis was conducted using Seurat clustering information and dimensionality reduction produced by UMAP. Cluster C0 was selected as start cluster, since C0 was mostly composed with normal chondrocytes, with specific parameters setting of “extend=‘n’, dist.method =“slingshot”. TradeSeq was used to discover transcription factors that were differentially expressed along trajectory lineages.¹⁹

Histology Analysis

Fresh human cartilage specimens from patients who underwent total knee arthroplasty were fixed in 4% paraformaldehyde overnight and decalcified in ethylenediaminetetraacetic acid (EDTA, G1105, Servicebio) for 3 weeks before embedding in paraffin. Then, 6 μm-thick sagittal paraffin sections were cut through the SIFK lesion, and the sections were stained with hematoxylin and eosin (H&E, Leagene), Safranin O/Fast Green (S2255-100 g, Sigma Aldrich) and Masson solutions (Servicebio) according to the manufacturer’s instructions. IHC staining was carried out using rabbit anti-TCF4 (1:50, C48H11, CST) as we described previously.²⁰

Statistical Analysis

The analysis was carried out in the R environment for statistical computing and visualization. For differential proportion, differential expression, and enrichment analyses, P values were adjusted for multiple hypothesis testing using the Benjamini–Hochberg method.

Results

Identification of Novel Chondrocyte Subpopulation in the SIFK Cartilage

Medial meniscus extrusion and insufficiency fracture of subchondral bone are believed to be the prime factors that contribute to abnormal stress loading on the medial cartilage and the development of SIFK.^8,21,22 Ultimately, this can lead to subchondral collapse and secondary OA. Until now, it has been unclear whether the transcriptional profiling in the cartilage of SIFK differ from those of OA. Herein, we sought to uncover the transcriptomic characteristics of cartilage obtained from SIFK using scRNA-seq, as illustrated in Figure 1A. The radiological changes and diagnosis were examined by weight-bearing X-ray and MRI (Figure 1B and Supplementary Figure S1). Of note, a typical subchondral bone lesion and medial meniscus extrusion were revealed by T2/T1-weighted MRI (Figure 1B, Koshino stage 2). Although conservative treatment has been suggested in the early stage of SIFK, in this case, UKA was chosen to rapidly alleviate the patient’s knee pain (Supplementary Figure S1). During the operation, cartilage specimens that covered the lesioned subchondral bone were dissected and subjected to chondrocyte isolation for scRNA-seq. The SIFK cartilage dataset was analyzed with a combination of normal and OA cartilage datasets deposited in GSE169454. After removing the batch effect by the Harmony package, the chondrocytes converged together but still exhibited a differential clustering pattern (Figure 1C and D, Supplementary Figures S2 and S3, Table S1). To further elaborate the cell subtype within the chondrocytes, we transferred the labeling of different stage OA chondrocytes identified by Ji et al¹⁴ onto the present integrated datasets. Several subpopulations, including effector chondrocytes (ECs), regulatory chondrocytes (RegCs), proliferative chondrocytes (ProCs) and hypertrophic chondrocytes (HTCs) were found in normal cartilage defined as C0, while ECs were hardly seen in both OA and SIFK group (Figure 1C and D). The C1 cluster was substantially composed of the OA cartilage with ProCs and prehypertrophic chondrocytes (preHTCs) as predicted cell subtype. Of significance, the overwhelming majority of C2 and C3 clusters were composed of the SIFK cartilage, C2 mainly contains the predicted homeostatic chondrocytes (HomCs), preHTCs and HTCs. C3 were predicted to be the fibrocartilage chondrocytes (FCs) and HTCs (Figure 1C and D). The differentially expressed genes (DEGs) in C0-3 clusters were further plotted and presented using a heatmap (Figure 1E). Hallmark, KEGG signaling pathway and GO term analyses of each cluster found that the activation of inflammatory signaling, including TNF-α signaling via NF-κB, IL-6/STAT3 signaling, and cytokine‒cytokine receptor interactions, was found in the OA cartilage, indicating a tendency of breaking down the extracellular matrix and collagens via chronic inflammation (Figure 1F and Supplementary Figure S4). However, inflammatory signaling was reduced in the SIFK cartilage; instead, robust activation of hallmarks of hypoxia and Wnt/β-catenin signaling as well as TGF-β signaling was found in clusters 2 (C2) and 3 (C3), respectively (Figure 1F). Furthermore, we used scMetabolism to investigate cartilage metabolism under different conditions. Unsupervised clustering revealed that normal, OA and SIFK cartilage cells had differential expression patterns of metabolism-associated genes and metabolic pathway scores (Supplementary Figure S5). C2 and C3 SIFK cartilage showed potential upregulation of glycolysis/gluconeogenesis and butanoate metabolism, which are highly relevant to hypoxia, sulfur and glycosaminoglycan metabolism, which are associated with chondrocyte differentiation or endochondral ossification activation. Collectively, these data suggest that a novel subpopulation of cells characterized by hypoxia and ossification activation was found in SIFK.

Figure 1 Single-cell RNA-sequencing (scRNA-seq) demonstrates reprogramming of SIFK cartilage. (A) Schematic of the experimental workflow used for scRNA-seq. (B) T2/T1-weighted MRI (coronal, sagittal, and axial views) showing a classic subchondral insufficiency fracture of the medial femoral condyle of the left knee joint. Red arrows show a diffuse high/low signal intensity in the medial femoral condyle. (C) Upper panel: UMAP plot of integrated datasets of chondrocytes from normal, OA, and SIFK cartilage, and the cell clusters are denoted as C0, C1, C2, and C3. Lower panel: Predicted identification from labeling transfer of GSE104782 (D) Upper panel: Split view of the integrated UMAP plot showing the differential clustering pattern of normal, OA, and SIFK cartilage. Lower panel: Split view of predicted identification from labeling transfer of GSE104782 (E) Heatmap illustrating the top 10 differentially expressed genes (DEGs) (p ≤ 0.05, fold change ≥ 1.5) in each cluster (C0-C3). (F) Heatmap showing the scaled pathway activities in each cluster. Abbreviations: Normal, healthy individual; OA, osteoarthritis patient; SIFK, Subchondral Insufficiency Fractures of the Knee.

Trajectile of Cartilage in SIFK Differs from That of OA

It is generally accepted that SIFK might lead to the late stage of OA. Therefore, it is presumed that SIFK cartilage may undergo a degenerative process similar to that of OA cartilage. Given the abovementioned novel subpopulation, we then performed trajectory analysis of these cells. Surprisingly, we observed that the trajectory of cartilage from normal conditions to OA was completely different from that of normal conditions to SIFK (Figure 2A and B). In addition, increased expression of biomarkers (COL2A1, ACAN, and SOX9) in hyaline cartilage was found in C2 and C3, along with elevated expression of fibrocartilage markers, including COL1A1, COL1A2, and COL3A1 (Figure 2C and D, Table S2). These findings are in accordance with previous histopathological observations that the formation of cartilage and fibrous tissue in the advanced stage of SIFK is indicative of delayed union or nonunion of the insufficiency fracture of the subchondral bone.⁴ Transcription factor analysis consistently revealed a shift to HIF1A-mediated hypoxia from normal to SIFK cartilage (Figure 2E–G). The Wnt/β-catenin pathway plays a key role in regulating chondrocyte hypertrophy, promoting bone formation and stimulating endochondral ossification, largely relying on the activity of LEF/TCF family members, including TCF4.²⁰ Herein, we found that TCF4 was predominantly expressed in the SIFK cartilage (clusters C2 and C3) (Figure 2E–G, Table S2). Furthermore, cartilage was obtained from OA (uninvolved lateral compartment of the knee as a normal control) and SIFK patients to validate the above results (Figure 2H, Supplementary Figure S6). As indicated, HE and S/O staining clearly revealed the reactive formation of cartilage and fibrous tissue and the activation of endochondral ossification, which was consistent with a previous report.²³ Accordingly, the expression of TCF4 was found to be increased in the ossification region by IHC staining (Figure 2H, Supplementary Figure S7). These results together suggested that transcription factors, including HIF1A and TCF4, might play a role in regulating hypoxia and endochondral ossification in SIFK.

Figure 2 Trajectory tracing of transcription factor alterations in the cartilage of SIFK. (A) UMAP plots revealing the lineage Slingshot trajectories of integrated normal, OA and SIFK cells; black lines represent predicted trajectories across pseudotime. (B) Slingshot trajectory plots of integrated normal and SIFK cells; red lines represent predicted trajectories across pseudotime (left). Lineage trajectories of normal and OA cell states; blue lines represent predicted trajectories across pseudotime (right). (C) UMAP plots showing chondrogenic marker gene (ACAN, COL2A1, SOX9) expression in normal, OA, and SIFK. (D) UMAP plots showing fibrosis marker gene (COL1A1, COL1A2, COL3A1) expression in normal, OA, and SIFK chondrocytes. (E) Pseudotemporal depiction of differentially expressed transcription factors starting from normal to SIFK (left) and OA (right). Values beyond the scale range are set to minimum or maximum. (F) Expression pattern of TCF4 and HIF1A along trajectories. Lineage 1: normal to SIFK; Lineage 2: normal to OA. (G) UMAP plots depicting the relative expression of HIF1A and TCF4 in normal, OA, and SIFK chondrocytes. (H) HE, Safranin O/Fast Green (S/O), and Masson’s trichrome staining of normal, OA, and SIFK tissue sections. TCF4 expression was examined by IHC. Red arrows indicate TCF4-stained cells. Scale bar: 200 μm.

Discussion

SIFK is characterized by acute onset of knee pain that mainly afflicts the elderly. It is often misdiagnosed in the clinic, especially in the early stage, since no change in appearance is observed on radiographs; however, lesions can become conspicuous as the disease progresses, manifesting as radiolucent areas with sclerotic halos and subtle flattening of the articular surface. In advanced cases, subchondral collapse and secondary OA can occur. Abnormalities of subchondral bone are strongly associated with cartilage degeneration in OA, which can be secondary to SIFK.⁵ Therefore, these progressive changes lead to a plausible assumption that cartilage in SIFK and OA undergo a similar degenerative process. In the present study, we reveal for the first time that cell type subpopulations and transcriptomic profiling of the articular cartilage in SIFK patients are completely different from that in normal and OA patients at the single-cell level.

Articular cartilage is a highly organized structure composed of multiple zones, namely, the superficial, middle, deep and calcified zones, which display three chondrocyte phenotypes: persistent, transient and hypertrophic cells. In recent studies, scRNA-seq has been used to identify several cell populations, including ECs, RegCs, HomC, ProCs, preHTCs, HTCs, FCs, senescent cells and cartilage progenitor cells, in OA.^11,14,24 In the current study, we further projected these nomenclatures onto normal and SIFK chondrocytes and shown that ECs which were known for association with metabolism reduced greatly in OA and SIFK chondrocytes. Of note, we were able to identify the overwhelming majority of C2 and C3 clusters in the SIFK cartilage, C2 were shown to contain HomCs and HTCs, while C3 were composed of FCs and HTCs, which is also completely different from OA and normal cartilage regarding hypoxia and endochondral ossification signaling (Figure 1). Our findings are in accordance with several previous studies demonstrating that the proliferative and redifferentiation abilities of chondrocytes from SIFK are largely different from those of normal and OA patients.^7–9

Substantial evidence has demonstrated the mutual regulation between cartilage and subchondral bone. For instance, during the epiphyseal extension process, chondrocyte-derived osteoprotegerin (OPG) facilitates bone formation by blocking osteoclast maturation through β-catenin signaling, and the senescence-associated secretory phenotype (SASP) of chondrocytes, which is characterized by the secretion of chemokines and proinflammatory cytokines, accelerates OA progression.²⁵ On the other hand, abnormalities in subchondral bone have been demonstrated to play a fundamental role in mediating cartilage degeneration.²⁶ In our study, the HE and S/O staining of SIFK cartilage revealed a reactive formation of cartilage and fibrous tissue and the activation of endochondral ossification (Figure 2), further validated the above scRNA-seq results (Figure 1). Mechanistically, chondrocyte trajectory analysis and IHC staining revealed that transcription factors, including TCF4, were highly upregulated during the occurrence of SIFK, which might drive the reactive formation of cartilage and fibrous tissue and the activation of endochondral ossification. This is supported by several studies demonstrating that TCF4 plays a key role in embryonic skeletal development and endochondral ossification.²⁷ Additionally, we found that several ossification-related genes, such as SOX5, TRPS1 and ZFHX4, were up-regulated in SIFK cartilage. SOX5 is reported to involve in inducing the expression of Sox9 targets and regulating endochondral ossification.^28,29 TRPS1 may coordinate the mineralization of the extracellular matrix by osteoblasts, thus regulating the development of bones formed by intramembranous or endochondral ossification.³⁰ Zfhx4 coordinates the transcriptional network of Osterix and thereby coordinates endochondral ossification.³¹ Therefore, it is likely that these transcription factors collectively promote the endochondral ossification during SIFK, which together supporting the concept of SIFK, rather than SONK.

Conclusion

To the best of our knowledge, this is the first report to elucidate the transcriptional changes in SIFK cartilage by using scRNA-seq. However, our study has several limitations. First, a limited number of samples were analyzed by scRNA-seq since this is a pilot study, which may lead to statistical bias caused by individual characteristics. Second, age- and stage-matched cartilage is required to eliminate any possible confounding factor. Therefore, if conditions permit, the number of samples from SIFK patients can be increased to confirm the precise mechanisms of chondrocyte metabolism in the SIFK population.