Establishment of a time-lapse chondrocyte dedifferentiation model

To establish a time-lapse chondrocyte dedifferentiation model, we cultured primary hyaline chondrocytes isolated from uncalcified cartilage in knee joints of newborn mice32,33 and replated each passage after doubling the cell quantity (seeding density: 2.4 × 104 cells per cm2; collecting density: 4.8 × 104 cells per cm2). The morphology of cells from passages (P) 0, 2, 4, and 8 was examined. During continuous culture, the chondrocyte morphology typically changed from round/polygonal to flattened/amoeboid-like, and the cellular size increased (Fig. S1A, B).

Col2-pd2EGFP reporter chondrocytes were used to examine the expression of the hyaline cartilage lineage marker Col2a1 during P0–8.34Col2a1 expression decreased as the passage number increased, which indicated a reduction in the functional phenotype (Fig. S1C, E). The transcriptomes of P0–8 chondrocytes were analyzed by performing scRNA-seq using a FluidigmTM C1 workstation (Fig. 1a). After the removal of low-quality data, we obtained the profiles of 634 cells, which covered 19 P0, 369 P2, 187 P4, and 59 P8 chondrocytes. A total of 6.9 ×107 reads and an average of 4 680 genes per cell were detected (GSE193742). The expression of representative dedifferentiation markers for P0–8 was detected and matched with the classical dedifferentiation patterns,8 including the reduced expression of cartilage functional markers (Col2a1, Acan and Col11a2) and the enhanced expression of fibrosis/degradation-related genes (Col1a1, Mmp13 and Mmp3) (Fig. S1F).

Fig. 1

ScRNA-seq identifies distinct cell subpopulations in the early and late passages. a The strategy for chondrocyte passaging and the scRNA-seq workflow. b T-distributed stochastic neighbor embedding (t-SNE) plotting of 4 chondrocyte clusters by Seurat unsupervised clustering; ProC proliferative cluster, EcmC extracellular matrix (ECM) cluster, MetC metabolic cluster, DegC degradative cluster. c t-SNE projection of cells and representative markers in 4 clusters. d Percentages of 4 clusters in different passages. e Heatmap of cluster-specific genes. f Relative expression level and percentage of representative markers in 4 clusters. g Enriched Gene Ontology (GO) terms in 4 clusters.

As a small number of cells was used to examine the transcriptomes of P0–8 chondrocytes, the results were compared with previously published scRNA-seq data for mouse primary chondrocytes and fibroblasts (GSE118236) (Fig. S1G).35 The principal component analysis (PCA) plot showed that all chondrocyte populations were distinctively distributed from mouse fibroblasts. Most P0 cells (18/19) merged with mouse primary chondrocytes from another data resource (GSE118236)35 (Fig. S1G), indicating that the evaluated P0 chondrocytes could represent typical primary chondrocyte subpopulations. We also confirmed the chondrogenic identity in the detected cells by demonstrating the high level of Col2a1/Col1a1 (ratio of chondrocyte marker expression to fibroblast marker) compared to that in mouse fibroblasts (Fig. S1H). This result was supported by the finding that the detection of chondrocyte dedifferentiation was not due to the overgrowth of fibroblasts present in the cultures.8 Overall, we established an in vitro cell model with typical phenotypic changes in chondrocyte dedifferentiation.

ScRNA-seq identifies distinct cell subpopulations in the early and late passages

Next, we analyzed chondrocyte subpopulations in P0–8 cells using Seurat unsupervised clustering from the R package.36 Four chondrocyte clusters were distinctly classified (Fig. 1b), with over 200 cluster-specific genes on average for each cluster (Fig. 1e, f, Table S1). We termed these clusters using Gene Ontology (GO) analysis (Fig. 1b). Cluster 1 was named the proliferative cluster (ProC, Fig. 1b) and was characterized by proliferative genes (Mki67 and Top2a, Fig. 1c, e, f). The ProC was enriched for GO terms that included the cell cycle and DNA replication (Fig. 1g). Cluster 2 was defined as the ECM cluster (EcmC), and it was characterized by highly expressed cartilage ECM genes, including Col9a1, Col11a2, and Col9a3 (Fig. 1c, e, f). The cluster was enriched for GO terms related to matrix production, including tRNA aminoacylation and proteinaceous ECM (Fig. 1g).

Cluster 3 was referred to as the metabolic cluster (MetC) and was characterized by metabolic-associated genes, retinol-binding protein 4 (Rbp4), isocitrate dehydrogenase 1 (Idh1), and secretory leukocyte peptidase inhibitor (Slpi) (Fig. 1c, e, f) and GO terms, including extracellular region, growth factor activity, and mitochondria. Cluster 4 was termed the degradative cluster (DegC) and contained metalloendopeptidase genes (Mmp3 and Mmp13, Fig. 1f). The GO terms extracellular matrix, actin cytoskeleton, and positive regulation of cell migration were enriched in the DegC (Fig. 1g).

To infer the dedifferentiation stages of each cluster, we plotted the cluster proportion of P0–8 (Fig. 1d). The ProC was mainly detected in P0 and P2 (42.1% of P0 and 26.9% of P2) (Fig. 1d). The EcmC had a similar distribution to the ProC (57.9% of P0 and 37.8% of P2) (Fig. 1d), which is likely to represent primary chondrocytes with abundant ECM features. The MetC was not detected in P0 and occupied less than 10% of P4 and P8; however, it accounted for 37.8% of P2 (Fig. 1d) and thus was assumed to be an early dedifferentiated cluster. The DegC was present in the majority of P4 and P8 (77.1% of P4 and 74.6% of P8) (Fig. 1d), which suggested that the DegC was a late-dedifferentiated cluster.

Pseudotime analysis reconstructs a loss-of-function trajectory

We reconstructed a pseudotemporal trajectory based on cell-to-cell similarity using the R package Monocle37 (Fig. 2a, b). Consistent with our previous hypothesis, EcmC cells were aligned to the beginning of the trajectory, while DegC cells were mostly aligned to the late stage. The unique MetC, which represented the early dedifferentiated stage, was predominantly postulated to be present in the midterm of the trajectory (Fig. 2a, b). In addition, we showed the patterns of the dedifferentiation and cartilage functional biomarkers (Mmp3, Mmp13, Adamts5, Col1a1, Col2a1 and Acan) on the trajectory axis as controls (Fig. 2c).

Fig. 2

Pseudotime analysis reconstructs a loss-of-function trajectory. a Pseudotime trajectory profiles of cells in the EcmC, MetC, and DegC revealed in a two-dimensional independent componentspace, reconstructed by Monocle. b Pseudotime trajectory profiles of cells in the EcmC, MetC, and DegC marked by passages. c Heatmap of representative dedifferentiation markers expressed in pseudotemporal order. d Expression pattern (left), enriched GOterms (middle) and examples (right) of dynamically regulated genes based on pseudotemporal order: immediately downregulated, gradually downregulated, delayed upregulated, gradually upregulated and tide wave style genes.

Next, we systemically characterized the time-lapsed transcriptome patterns by listing the top 1 000 genes that were dynamically regulated (q < 0.01). The genes were classified into five groups based on their expression patterns (Fig. 2d, Table S2): immediate downregulation, gradual downregulation, delayed upregulation, gradual upregulation, and tide wave pattern. Briefly, genes that were expressed early and showed immediately downregulated expression were assigned to GO terms that included proteinaceous ECM, tRNA aminoacylation, and ribosome (Fig. 2d). Genes that showed gradually downregulated expression were aligned with the GO terms translation, small ribosome subunit, methylosome, and protein import into the mitochondrial matrix (Fig. 2d). In contrast, genes with delayed upregulated expression were involved in the inflammatory responses cellular responses to calcium (Ca2+), cellular responses to interleukin (IL)-1, and cytokine activity (Fig. 2d). Genes involved in collagen fibril organization, collagen trimer, blood vessel development, and metalloendopeptidase activity showed gradually upregulated expression during chondrocyte dedifferentiation (Fig. 2d). Transcription factor (TF)-coding genes with upregulated or downregulated expression were also detected (Table S3), and those regulating transcription, Ca2+responses, skeleton morphologies, and cAMP responses, including Jun (C-JUN, AP-1 transcription factor subunit), Fos, Runx1, and Creb3l1 (c-AMP responsive element binding protein 3 like 1), were highly expressed at the late stages (Fig. S2A, B). The TF genes with downregulated expression were involved in DNA binding and DNA stability38 and were associated with the minichromosome maintenance complex (Fig. S2A, C).

In all, genes responsible for ECM production and protein translation showed widely downregulated expression during chondrocyte dedifferentiation, while the genes for stress, inflammatory responses, and collagen reorganization showed significantly upregulated expression.

Early dedifferentiated chondrocytes highly express a wide spectrum of metabolic genes

We further characterized early dedifferentiated chondrocytes. The cluster distribution indicated that P2 cells partially retained the features of P0 chondrocytes (ProC and EcmC) (Fig. 1d) but also contained the unique subpopulation MetC, which mostly represented early dedifferentiated chondrocytes (Fig. 2a). The MetC was enriched for metabolism-associated genes, including Rbp4 (encoding adipokine retinol-binding protein 4), Idh1 (encoding isocitrate dehydrogenase 1) and Mpv17 (encoding mitochondrial inner membrane protein 17) (Fig. 3a).

Fig. 3

Early dedifferentiated chondrocytes highly express metabolic genes and obtain a glycolytic phenotype. a Typical cluster marker expression based on the pseudotemporal order. b Representative genes of glycolytic processes and cellular responses to oxidative stress in the EcmC, MetC, and DegC. c Oxygen consumption rate (OCR) assays of P0-8 chondrocytes. d Extracellular acidification rate (ECAR) analysis of P0-8 chondrocytes. e Two-dimensional bioenergetic profiles of P0-8 chondrocytes. f, g Other test parameters of OCR and ECAR assays of P0-8 chondrocytes. h Relative total intracellular ROS production in P0-8 chondrocytes. i, j Relative mitochondrial ROS production in P0-8 chondrocytes and the quantitative data, detected by MitoSOX Red staining; scale bars: 50 μm. All data are the mean ± SEM. *P < 0.05, **P < 0.01, ***P < 0.001.

The GO terms enriched in the MetC relative to the EcmC included extracellular region, sterol biosynthetic process, chemotaxis, glutathione transferase activity, and mitochondrial outer membrane (Fig. S3B). KEGG pathway analysis suggested that genes assigned to the cytochrome P450 metabolism and glutathione metabolism pathways had significantly upregulated expression in the MetC (Fig. S3C), with the representative genes being Mgst1, Mgst3 (genes encoding microsomal glutathione S-transferase 1 and 3), Gstt1, Gstt2 (genes encoding glutathione S-transferase theta 1 and 2), Cyp51 (cytochrome P450, family 51), and Sod3 (superoxide dismutase 3) (Fig. S3D). Furthermore, we revealed that genes involved in the glycolytic process (Cst3, Prdx5, Sod3, Gpx1 and Pebp1, etc.) and cellular response to oxidative stress (Aldoa, Gpi1, Pgk1, Gapdh and Eno1, etc.) were activated in the MeC compared to the EcmC and DegC (Fig. 3b). Thus, the major transcriptional features of the MetC have been demonstrated to be involved in metabolic processes.

In addition, genes associated with chemotaxis and chondrocyte hypertrophy were highly expressed in the MetC (Fig. S3E, G). Consistently, approximately 39% of P2 chondrocytes were shown to have a migratory phenotype with a fibroblast-like shape in comparison with the typical round shape predominantly present in primary chondrocytes (Figs. S1A and S3H). Genes associated with proliferation did not show significantly upregulated expression in the MetC (Fig. S3F). Thus, the results suggested that the genes typically expressed in the early dedifferentiated cluster are associated with metabolic activation, chemotaxis, and chondrocyte hypertrophy but not proliferation.

Early dedifferentiated chondrocytes show a glycolytic phenotype during metabolic reprogramming

The unique subpopulation MetC, characterized by metabolic activation, occupied distinct proportions in different passages (Fig. 1d) and likely represented early differentiated chondrocytes. Thus, we further examined whether metabolic activation could be validated with cellular metabolic behavior. We performed live-cell metabolic assays to measure the oxygen consumption rate (OCR) and extracellular acidification rate (ECAR) in P0–8 chondrocytes (Fig. 3c–g). Briefly, the results indicated an anaerobic-to-aerobic metabolism switch during chondrocyte dedifferentiation (Fig. 3c–g). P0 and P2 cells had a significantly higher glycolytic rate and capacity than P4 and P8 cells (Fig. 3d, g, P < 0.001), while P4 and P8 cells showed increased maximal OCR respiration after FCCP injection (P < 0.001), suggesting an enhanced oxidative phosphorylation capacity to meet the increased energy demand (Fig. 3c, f). Increased levels of intracellular total reactive oxygen species (ROS) were also detected in later passages, especially in P8 chondrocytes (Fig. 3h and Fig. S3J, P < 0.001). The analysis specifically measuring mitochondrial ROS showed that late-stage chondrocytes had an increased level of mitochondrial ROS, in a similar trend to that of total ROS (Fig. 3i, j).

Notably, P2 cells showed a glycolytic phenotype when compared to cells at other passages, according to the metabolic phenotype map (Fig. 3e), and had a higher level of basal ECAR and glycolysis rate than P0 cells (Fig. 3g). The OCR results demonstrated that P2 chondrocytes produced less mitochondrial ATP than the cells from other passages, with more proton (H+) leakage uncoupled to ATP production (Fig. 3f). Increased H+ leakage has been documented as a protective mechanism against oxidative stress.39,40 Consistent with this finding, ROS production did not significantly increase in P2 chondrocytes compared to P0 chondrocytes (Fig. 3h, Fig. S3J, I, J).

Overall, live-cell metabolic assays demonstrated that dedifferentiated chondrocytes underwent a metabolic transition with a higher aerobic respiration capacity. However, P2 cells showed a unique glycolytic metabolic phenotype with increased mitochondrial proton uncoupling and reduced ROS than the cells at other passages. This finding was consistent with the result that P2 contained the highest population of the MetC, a cluster with activated genes for the glycolytic process and cellular response to oxidative stress.

Late dedifferentiated chondrocytes are characterized by ultrastructural changes of metabolic stress

Next, we focused on late-dedifferentiated chondrocytes. The DegC represented the majority of P4 and P8 chondrocytes (Fig. S4A, B) and contained genes involved in cytoskeleton reorganization, matrix degradation, inflammatory response, and chondrocyte hypertrophy (Fig. 1g, Fig. S3G, I). In particular, the DegC was also characterized by features of mitochondrial stress, as well as the mitochondrial unfolded protein response,41,42 a self-protective mechanism that promotes cell survival against external stimuli and delays mitochondrial dysfunction.41,43,44,45 Relevant characteristics for the mitochondrial unfolded protein response, including attenuated translation, impaired mitochondrial homeostasis, activated mitochondrion-to-nucleus retrograde pathways, and epigenetic alterations, have been reported.43,44,46 Consistently, in our results, the number of genes expressed per cell was reduced in the DegC (Fig. S4C), consistent with a previous study,18 and the lowest expression of translation-relevant genes was detected (Fig. S4D). Genes related to mitochondrial function (translocase and porins), including Timm17a, Timm23, Tomm22, Vdac2 and Vdac3, showed slightly upregulated expression in the MetC but downregulated expression in the DegC (Fig. S4D). Genes of the MAPK pathway and AP-1 family that are involved in mitochondrial retrograde signaling47,48 and are known nuclear responsive factors were highly expressed in the DegC (Fig. S4D).

To confirm impaired mitochondrial homeostasis in late dedifferentiated chondrocytes, we examined the ultrastructure of P0-8 chondrocytes by transmission electron microscopy (TEM) (Fig. 4a). Quantitative analysis demonstrated an increased number of mitochondria in P4 and P8 cells (Fig. 4a, d), matching their aerobic phenotype. Specifically, P0 displayed more homogeneous and round mitochondria, with lower aspect ratios (Fig. 4a, b, d), while P2 mitochondria were relatively elongated (Fig. 4a, b, d) compared to those of other groups. In contrast, P4 and P8 contained perinuclear, heterogeneous mitochondria with disorganized configuration and an electron-dense damaged structure49 (red arrows) (Fig. 4a, b, d). Similarly, MitoTracker images showed swelling and fragmentation of mitochondria in late dedifferentiated cells (Fig. 4c).

Fig. 4

Late dedifferentiated chondrocytes exhibit ultrastructural changes in metabolic stress. a Representative images of mitochondria in P0-8 chondrocytes by transmission electron microscopy; scale bars: 1 μm; blue arrows: mitochondria, red arrows: damaged mitochondria, green arrows: endoplasmic reticulum. b Schematic diagram of mitochondrial morphological features in primary (P0), early dedifferentiated (P2) and late dedifferentiated chondrocytes (P4 and P8). c Representative images of MitoTracker-stained mitochondria in P0-8 chondrocytes; scale bars: 20 μm. d Quantitative analysis of mitochondrial TEM images of P0-8 chondrocytes. e, f Representative images and quantitative analysis of DAPI-stained nuclei in P0-8 chondrocytes by confocal microscopy; scale bars: 20 μm. g Representative genes with P8-specific open regions in promoters (promotor regions are marked in green). h Representative enriched TF binding motifs in P8. i Heatmap of enriched motifs in P0 and P8. All data are the mean ± SEM. *P < 0.05, **P < 0.01, ***P < 0.001.

Overall, the results showed that late dedifferentiated chondrocytes contained an increasing quantity of mitochondria to meet metabolic demand. However, structural damage to the mitochondria was observed. As mitochondrial morphology is controlled by metabolism,50 this finding may indicate the dynamic fusion and fission of mitochondria.

Late dedifferentiated chondrocytes exhibit stress-associated chromatin remodeling

Together with the organelle structure, nuclear morphology changed dramatically with increased passage number. High-resolution confocal images showed that nuclei in P4 and P8 chondrocytes were significantly enlarged with enhanced integrated fluorescence intensity (Fig. 4e, f), suggesting a structural alteration in nuclear architecture and chromatin condensation.51

Next, an assay for transposase-accessible chromatin with high throughput sequencing (ATAC-seq) was performed based on a previous method52 (P0 vs. P8, 2 replicates for each, Fig. S4E) to confirm the chromatin landscape modifications in late-dedifferentiated chondrocytes. A decrease in global chromatin accessibility was detected in P8 chondrocytes (Fig. S4F, GSE193743), and stage-specific peaks at 16754 and 5971 were identified at P0 and P8, respectively (Fig. S4G). P8 cells were shown to have fewer open regions near gene promoters (<1 kb) and fewer TF binding loci near transcriptional start sites (TSS, <1 kb) (Fig. S4I–L). In contrast, the proportion of distal intergenic regions in P8 cells decreased from 43.3% to 37.5% in P0 cells (Fig. S4K).

Interestingly, similar epigenetic alterations have been reported in stress responses of other lineages.53,54 The transcription of most cellular genes was disrupted, but specific downstream effectors were activated. Specifically, we identified representative P8-specific regions with enhanced accessibility near promoters and found that the genes associated with these regions were relevant to fibrosis (Col1a1, S100a4, and Fbln2), immune responses (Ccl2, Ccl2 and Ly6c1), and metabolic stress (Gstm2 and Sod3) (Fig. 4g). We also compared the TF binding motifs in P0 and P8 cells (Fig. 4h, i). The data revealed that binding sites for TFs, which play essential roles in mitochondrion-to-nucleus signaling, oxidative stress, and inflammatory responses, including the ATF3, BATF, and AP-1 family, were strongly enriched in P8 cells (Fig. 4h, i).

Thus, late dedifferentiated chondrocytes exhibited epigenetic remodeling with a global decrease in chromatin accessibility and increased open regions in genes responsive to metabolic stress.

Manipulating mitochondrial F1Fo-ATPase efficiently ameliorates early dedifferentiation

Having identified the increasing metabolic stress and mitochondrial dysfunction during chondrocyte dedifferentiation, we hypothesized that manipulating stress signaling might influence the chondrocyte functional phenotype. F1Fo-ATPase is the prime ATP producer on the inner mitochondrial membrane and essentially regulates mitochondrial function.55 This molecule synthesizes ATP by pumping H+ across the inner mitochondrial membrane to the mitochondrial matrix. In some pathological events, when mitochondrial homeostasis is compromised, F1Fo-ATPase can degrade ATP and reverse the direction to generate proton backflow (from the mitochondrial matrix to the intermembrane space).56,57,58 The chemical compound BTB06584 (BTB) selectively inhibits F1Fo-ATPase ATP hydrolysis and reduces F1Fo-ATPase-driven H+ backflow without compromising ATP synthesis58 (Fig. 5a).

Fig. 5

Manipulating mitochondrial F1F0ATPase efficiently ameliorates early dedifferentiation. a Schematic diagram of BTB06584 (BTB)’s effect on F1Fo-ATPase. b Representative immunostaining images of COL2 and MMP13 in the BTB-untreated (negative control, NC) and treated P2 chondrocytes; scale bars: 50 μm. c, d qPCR detection of representative dedifferentiation genes in the NC and BTB-treated P2 chondrocytes. e Quantitative data of COL2 and MMP13 immunostaining of the NC and BTB-treated P2 chondrocytes. f Quantitative analysis of the cellular and nuclear sizes of the NC and BTB-treated P2 chondrocytes. g Quantitative analysis of relative ROS production in the NC and BTB-treated P2 chondrocytes. h Quantitative analysis of relative mitochondrial ROS production in the NC and BTB-treated P2 chondrocytes. i Heatmap of differentially expressed genes and enriched GO terms in the BTB-treated vs. NC chondrocytes, detected by bulk RNA-seq. j Up- and downregulated KEGG pathways in the BTB-treated P2 vs. NC chondrocytes. The final concentration of BTB was 2.5 μmol·L−1. “NC” (negative control) in this figure represents the treatment with 0.025% dimethyl sulfoxide (DMSO), the solvent of BTB. All data are the mean ± SEM. *P < 0.05, **P < 0.01, ***P < 0.001.

BTB treatment of P2 and P4 chondrocytes resulted in different phenotypic changes. Immunostained images demonstrated that BTB significantly promoted COL2 expression and inhibited MMP13 expression in P2 chondrocytes (Fig. 5b, e). We then conducted real-time qPCR, and the results showed higher expression levels of representative functional genes (Acan, Sox6, Sox9 and Col2a1) and suppression of matrix degradation genes (Adamts5 and Mmp13) in the BTB-treated group (Fig. 5c, d). The BTB-treated P2 chondrocytes did not exhibit distinct cellular or nuclear sizes (Fig. 5f). We also compared the phenotypes of the BTB-treated P2 cells with those of the BTB-untreated P1 and P2 chondrocytes. The immunostaining results of COL2, ACAN and MMP13 demonstrated that BTB treatment could not fully revert P2 to P1 but rescued most of the phenotype loss (Fig. S5A, B).

In live-cell metabolic assays, P2 cells showed a high level of H+ leakage into the mitochondrial matrix, which was assumed to be due to ROS production,39 one of the triggers of mitochondrial-mediated stress.59,60,61,62 As BTB inhibited ATP hydrolysis and thereby reduced H+ backflow, we hypothesized that BTB's effect might be associated with decreased ROS production in P2 chondrocytes. We measured both total intracellular and mitochondrial ROS (Fig. 5g, h, Fig. S5D, E). Consistently, the results showed that the BTB-treated P2 chondrocytes significantly produced less ROS. No morphological difference was observed in the mitochondria of the BTB-treated and control P2 cells (Fig. S5C).

In P4 cells representing late dedifferentiated chondrocytes, BTB significantly inhibited MMP13 expression but had no effect on COL2 expression (Fig. S6A, B). BTB significantly inhibited only the matrix degradation genes ADAMTS5 and MMP13 (Fig. S6C, D). The BTB-treated P4 chondrocytes had a relatively smaller cellular size (Fig. S6G) but a similar nuclear size to the control group (Fig. S6G). Moreover, in P4 chondrocytes, ROS production was not suppressed by BTB (Fig. S6E, F), and the fragmentation and swelling of mitochondria was not significantly altered (Fig. S6H). TEM revealed that damaged and heterogeneous mitochondria were not obviously reduced in the BTB-treated P4 cells (Fig. S6I).

Thus, we revealed that the chemical inhibitor BTB could effectively ameliorate functional phenotype loss in early dedifferentiated chondrocytes (P2) but not in late dedifferentiated chondrocytes (P4). ROS production could be significantly suppressed in P2 chondrocytes. The mitochondrial structural damage in P4 cells was not obviously ameliorated.

BTB ameliorates early phenot