Two reactive behaviors of chondrocytes in an IL-1β-induced inflammatory environment revealed by the single-cell RNA sequencing

Identification of inflammation-responding cell clusters and their expression signature

The intact cartilage was isolated and digested from discarded finger tissue after excision surgery from the same donor. The in vitro expanded chondrocytes were then subjected to designed experiments. Figure 1A shows a schematic showing the process of species access and data collection along with which cell proliferation and morphology were also explored. Enhanced catabolism and induced inflammatory in chondrocyte under IL1β stimulation were confirmed by qPCR which showed up regulated gene expression of MMP1, MMP3, MMP13 and ADAMTS5, and genes of inflammatory factors like TNF and IL6 (Supplementary Figure 1). The well-distributed cytoskeleton was observed, and several cellular pseudopodia could be identified, although they were not quite as distinct, in both articular chondrocytes (ACs) treated with IL-1β and controls (Supplementary Figure 2A). Correspondingly, no significant difference in proliferation activity was identified in ACs treated with IL-1β and control ACs, after a 3-day culture (Supplementary Figure 2B). During a short-term inflammatory stimulation period, ACs demonstrated few changes in proliferation rate or macroscopic structure. Thus, the following analysis would mainly focus on the function transformation of ACs.

Figure 1. Work flow and identification of cell clusters. (A) A schematic of the process of species access and data collection. (B) The u-map analysis of the integrated data resulted in the identification of nine distinct cell clusters. (C) The distributions of the cell clusters on each sample; different colors represent different cell clusters and each plot indicates one cell. (D) The stacked chart of the proportions of the distinct cell clusters in the three samples. (E) The distribution ratio and quantity in samples of each cell cluster. The colors of the outer circle of each pie chart represent cell clusters. Each pie chart illustrates the ratio of one certain cluster in three samples and the accompanying histogram shows the quantity change of one certain cluster in three samples. (F) The graphical heatmap of the top 5 differentially expressed genes for each cell cluster.

After quality control, 4,703, 6,170, and 7,348 cells were recruited in downstream analysis after individually sequenced using scRNA-seq to characterize their transcriptomes. The u-map analysis of the integrated data resulted in the identification of nine distinct cell clusters, with numeric distinction among separate subpopulations (Figure 1B). The separation distance between all the identified cell clusters could easily distinguish the control cells and IL-1β stimulated cells among which the former was mainly consisted of clusters 2 and 3 and the latter of clusters 0, 1, 4, 5, 6,7. This was further confirmed by the split plot of cells from each sample (Figure 1C) which showed that clusters 1, 5 and 7 enriched in the 24-h sample, and clusters 0, 4, and 6 were mainly from the 48-h sample while the cluster 8 consistently existed in all the samples. Here, we defined the cell cluster by different gene expression patterns in order to explore the function transformation from the initial to the terminal phase.

The two dominating clusters in control sample shifted dramatically into the six clusters in stimulated samples with significant increases or decrease in the proportions of certain clusters following the IL-1β induction of inflammation in the culture environment (Figure 1D, 1E). Specifically, the proportions of cells grouped in clusters 1,7, and 5 increased sharply, from baseline at 0-h to highest quantity level in the 24-h sample, whereas the proportions of cells grouped in clusters 2,3 almost disappeared after IL-1 stimulation. Subsequently, clusters 0,4, and 6 gradually took place of the previously dominating clusters 1, 5 and 7 along with the continuing inflammation environment and constituted the main composition in the 48-h sample. So, the clusters 0,4, and 6 may serve as the potential inflammation-special cell cluster in the terminal phase.

There is a continuous but heterogeneous dynamic changes of cell function among these cell clusters along with the IL1β stimulation. Differential gene expression analysis was conducted among the u-map-separated clusters, and the graphical representation (heatmap) indicated the top 5 differentially expressed genes for each subpopulation (Figure 1F). Notably, although the gene expression patterns were distinct among clusters with marked fluctuations in cell quantities, the specifically high-expressed genes among clusters 0, 4 and 6 were most commonly associated with the metabolism of the ECM such as MMP3, COL1A1, COL1A2, and COL3A1 and inflammatory factors like IFI27, TNFAIP6, TGFBI. Interestingly, clusters 1,5 and 7 shared some common high-expressed genes with clusters 0,4, and 6, most of which were relating to collagen synthesis. Although these cell clusters exhibited dynamic changes in terms of cell quantities and genetic functions in different samples, cluster 8 showed a relatively stable status independent of stimulation whether in cell numbers or gene expression pattern.

The correlation between clusters-special genes and functional gene sets and previously-identified marker genes of subpopulation of chondrocytes

To classify the function of each cluster, we first referred pathway enrichment based on gene set variation analysis, and then calculated the correlations between the expression signature of the 9 clusters and constructed functional gene sets. Pathway enrichment analysis revealed an increasing relation to inflammatory response from clusters in original state like cluster 2 and 3 to the terminal clusters such as clusters 0 and 4 while clusters 1,7, and 0 closely linked (Figure 2A). The heatmap of correlation between cluster-specific expressed genes and constructed functional gene sets also indicated that cluster 1-specific genes were highly related to the inflammatory gene set among the major cell clusters in stimulated samples. However, gene sets containing the genes that regulate circadian rhythm, proliferation, or microvasculature were not significantly associated with any cluster-special genes (Figure 2B).

Figure 2. Correlation between expressions of clusters-special genes and functional gene sets and previously-identified marker genes of subpopulation of chondrocytes. (A) The heatmap of enriched pathways in clusters by GSVA. (B) The bubble chart of correlations between clusters-special genes and previously-identified marker genes of subpopulation of chondrocytes. (C) The heatmap of correlations between subpopulation of chondrocytes and constructed gene sets.

The subpopulation of chondrocytes was previously identified by as hypertrophic chondrocytes (HTCs), prehypertrophic chondrocytes (preHTC), fibrocartilage chondrocytes (FC), regulatory chondrocytes (RegC), homeostatic chondrocytes (HomC), effector chondrocytes (ECs), and proliferative chondrocytes (ProCs). Correlation between cluster-specifical genes expression and previously-identified marker genes expression of subpopulation of chondrocytes showed that cluster 1 was associated with preHTC and FC (cor = 0.82), cluster2 was highest related with RegC (cor = 0.83) and cluster 3 with HomC (cor = 0.60) and RegC (cor = 0.68), cluster4 with HomC(cor = 0.60) and RegC (cor = 0.63), cluster5 with RegC (cor = 0.77) and HomC (cor = 0.64), cluster 8 with preHTC (cor = 0.66) (Figure 2C and Supplementary Figure 3). Notably, earlier appeared clusters 2,3 and 5 were all has a similarity with RegC and HomC while the terminal clusters 0 and 7 have no obvious relationship with any predefined subpopulations. The specific correlation heatmap was presented in Supplementary Figure 3.

The response-to-inflammation paths constructed by single-cell trajectory analysis

A trajectory, with three main branches and one connecting points, was constructed by Monocle 2, in which cells were grouped into three states (Figure 3A and Supplementary Figure 4A). The trajectory resulted in two paths along with that cells transformed differently, which was further validated below. We found the cells in clusters 2 and 3 to be distributed in the earlier stage of the trajectory, which covered state 1 and was enriched with most of the cells from control sample, whereas the cells in clusters 1,7, and 0 were located in the later phase and extend into the ending along one branch (Supplementary Figure 4B and Figure 3B). Notably, the upregulated genes in cell clusters 1 and 7 consistently enriched in pathways and functional gene sets correlated with inflammation and ECM as depicted above (Figure 2A, 2B) and were primarily enriched in states 2 corresponding with one of the end directions of the pseudo-time trajectory. Thus, the increasing enrichment of end-state cells, the cluster-0 cells, were considered to be the final inflammation-responded chondrocytes during the late inflammatory stage. Therefore, we termed the path, in which cells were distributed from the root state 1 to the end state 2, the “response-to-inflammation path.”

Figure 3. Cell distribution along pseudo-time trajectory and construction of the response-to-inflammation path. (A) A trajectory with three branches and one connecting points was constructed by Monocle 2. (B) The time-line of the pseudo-time trajectory; the darker blue indicates earlier pseudo-time. (C) The distribution of cluster2, 1, 7, and 0 along the pseudo-time trajectory. These cells compromised a complete evolutionary route without branches. (D) The distribution of cluster3, 5, 6, and 4 along the pseudo-time trajectory. These cells compromised a different evolutionary route from the former with differentiated branches. (E–G) The expression level of specific representative genes of terminal cluster 0 in cluster2, 1, 7, and 0 are highlighted along the response-to-inflammation path and real-time line. (F–H) The expression level of specific representative genes of terminal cluster 4 in cluster3, 5, 6, and 4 along the atypical response-to-inflammation path and real-time line.

When examining changes in the transcript abundance of the inflammation genes in each cluster along the peseudotime trajectory, gene expression was observed to fluctuate with IL-1β stimulation but the trend was observed in one path. Specifically, the gradually increasing expression of BST2, IFI27, ISG15, MMP3, and TNFAIP6, most of which were inflammatory factors, were not only reflected by the order of clusters 2,1, 7, and 0, namely the “response-to-inflammation path”, but also match the time line of stimulation (Figure 3E, 3G). Namely, this peseudotime trajectory was consistent with the actual time line designed in this study.

The other path illustrated the evolving route from root cells of cluster 3 to the other ending cells of clusters 5, 6, and 4, of which the enriched pathways have relatively little relationship with inflammation and ECM metabolism (Figure 2A, 2B). Considering these cells ended at a different state comparing with the above “response-to-inflammation path”, this trajectory was termed as “atypical response-to-inflammation path”. Correspondingly, highly expressed genes of the terminal cluster 4 such as CDK1, HIST1H1B, MKI67, MX1, and TOP2A, showed a significant variation but not an exact trend along the “atypical response-to-inflammation path” both in pseudo-time and real-time (Figure 3F, 3H).

The characteristics of response-to-inflammation paths validated by gene function annotation and the identification of the key genes depending cell fate

To explore changes in cell function as cells evolve into the final state along this trajectory, we evaluated the gene expression pattern associated with point 1 where cells separately begin moving into different states, one of which mainly contained the inflammation amplification clusters and the other was atypical response-to-inflammation path. The expressed genes associated with cell fate at point 1 were shown in Figure 4A. Specifically, cells with up-regulated genes such as CHI3L1, CLU, IGFBP5, CCL20, G0S2, MMP13, MMP1, BST2,C15orf48, COMP, SAT1, TNFAIP6, SERPINE2, LUM, and MMP3 were assigned to cell fate1 which also represented the direction of the response-to-inflammation path while those with up-regulated genes like CCNF, FAM83D, HIST1H1D,BRCA2, YWHAH, DNMT1, ANP32E, and CENPH were allocated to cell fate2 following the atypical response-to-inflammation path. Consistently, the fate1 cells were mainly constructed by clusters 0, 1, and 7 while the fate2 cells mostly obtained clusters 4, 5, and 6 (Figure 4B, 4C), confirming the existence of the defined two paths through which chondrocytes reacted differently (Figure 4D).

Figure 4. Identification and analysis of genes involve in the connected point in the pseudo-time trajectory (A) The expression patterns of the genes associated with cell fate at point 1. Cell fates depend on the up- or down-regulation of these genes. (B, C) The pie chart of the proportions of the distinct cell clusters in cell fate 1 and 2. (D) The exhibition of cell transformation along two defined pathways. The grey arrows indicate the response-to-inflammation path and the azure arrows indicate the atypical response-to-inflammation path. (E, F) Significant signaling pathways and biological processes identified by Gene set enrichment analysis (GSEA) and gene ontology (GO) analysis conducted on fate 1 (the response-to-inflammation path) associated genes. (G, H) Significant signaling pathways and biological processes identified by Gene set enrichment analysis (GSEA) and gene ontology (GO) analysis conducted on fate 2 (the atypical response-to-inflammation path) associated genes.

The top pathways enriched in fate1 associated genes showed a corresponding results and are listed in figure 4E–4H. Notably, significantly enriched pathways relating to complement and allograft rejection were suggested that immune system may also play a role in the cartilage degeneration apart from the common inflammation regulatory pathways, such as the tumor necrosis factor (TNF) signaling pathway via NFKB, up-regulated KRAS signaling and the IL2/STAT5 signaling pathway. And pathways relating to apoptosis and reactive oxygen species indicated more specific biological processes involved in the cell’s transformation along the response-to-inflammation path. However, top pathways enriched by fate 2 associated genes mainly participated the cell cycle regulation such as E2F target, G2M checkpoint and mitotic spindle, representing a completely different cell state from cells in fate 1.

A gene ontology (GO) analysis conducted on these fate-depending genes is depicted in Figure 4F, 4G. Specifically, genes regulating cell transformation toward fate1 were more responsive to extracellular modification, with the most significantly enriched genes associated with extracellular structure organization during biological processes, and other classical phenotype in OA chondrocytes such as the ossification, connective tissue development, cartilage development, and collagen fibril organization, all involved in ECM catabolic functions. In comparison, genes regulating cells to the state 3 were more enriched in cell proliferation like nuclear division, organelle fission, chromosome segregation and DNA replication.

Ethics statement

Investigation has been conducted in accordance with the ethical standards and according to the Declaration of Helsinki and according to national and international guidelines and has been approved by the Committee for Medical Research Ethics of Tongji Hospital.

Isolation and culture of human AC

ACs, obtained from an interphalangeal joint of a polydactyl child, were isolated from discarded finger tissue after excision surgery. All cells used during single-cell sequencing were obtained from the same donor. The isolated cartilage was intact, containing all cartilage zones, and was used within 6 h of dissection. Briefly, the cartilage was shredded into small pieces, digested with collagenase type II (Sigma-Aldrich, St. Louis, MO, USA) for no more than 4 h, and seeded in culture medium consisting of DMEM/F12 (Thermo Fisher Scientific/ Gibco, Waltham, MA, USA) and containing 10% human platelet lysate-rich plasma (hPLP), 100 units/mL penicillin, and 100 μg/mL streptomycin (Sigma-Aldrich). ACs were cultured for 13 days (three passages). The obtained cells were reseeded in six-well plates and treated with human IL-1β (h-IL-1β) (10 μg/mL) for 0, 24, and 48 h. Informed consent has been obtained. All experiments were conducted under aseptic conditions.

Cell proliferation and morphology

Cell counting kit 8 (CCK8, Boster Biological Technology Co., CA, USA) was used to test the proliferative activity of ACs cocultured with h-IL-1β (10 μg/mL) for 24, 48, and 72 h. Approximately 2 × 10⁴ ACs, in 1 mL DMEM/F12, were seeded in a 24-well plate and cultured at 37° C, in a 5% CO₂ environment, for 24, 48, and 72 h. The culture medium was refreshed every other day. The absorbance of the final reaction solution, which was prepared according to the manufacturer’s instructions, was measured using a spectrophotometer (Bio Tek Instruments Inc., VT, USA).

Phalloidin (Sigma-Aldrich, USA) and 4′,6-diamidino-2-phenylindole (DAPI), Boster Biological Technology Co., CA, USA) were used to the visualize cell skeleton and cell nuclei, respectively. Briefly, ACs in 24-well plates treated with h-IL-1β (10 μg/mL) for 0, 24, and 48 h rinsed with phosphate-buffered saline (PBS), and fixed with 4% paraformaldehyde. Then, they were consecutively stained with phalloidin for 30 min and DAPI for 15 min in the dark. A fluorescence microscope (Nikon Ti-S, Tokyo, Japan) was used to observe changes in cell morphology, after air drying at room temperature.

The catabolism-related and inflammatory gene expression

Four critical catabolism-related genes, MMP1, MMP3, MMP13 and ADAMTS5, and two inflammatory genes, IL6 and TNF, were examined. The quantitative reverse transcription PCR (qRT-PCR) were applied to measure the gene expression levels. Briefly, cells were cultured in 12-well plates treated with h-IL-1β (10 μg/mL) for 0, 12, 24, and 48 h. The total cell biomass from three randomly selected samples in each timepoints were lysed with Trizol (Invitrogen, Thermo Fisher, MA, USA), and the Superscript II first-strand cDNA synthesis kit (Takara, Beijing, China) were used to reverse transcribe RNA. The SYBR Green (TOYOBO Biotech Co., Osaka, Japan) kit was used to prepare PCR mix and a Bio-radMyiQ2 thermocycler (Bio-Rad Laboratories Inc., Singapore) was used for quantification. All the primers that were used are listed in Table 1.

Cell collection and 10x genomic sequencing library construction

The obtained cells reseeded in six-well plates and treated with h-IL1β (10 μg/mL) for 0, 24, and 48 h were harvested at each timepoint. The collected cells were centrifuged at a speed not exceeding 300 × g and then washed three times in PBS. The cell precipitation was resuspended by 1 mL 1x PBS containing 0.04% bovine serum albumin and filtered through a 40-μm cell strainer. After washing twice, 100 μL PBS was added to the cell precipitate, to obtain a single-cell dispersion suspension containing a sufficient number of cells. Cells treated with the Single-Cell 3′ Reagent Kit User Guide v2 (10x Genomics) were then transferred to a Chromium Controller instrument (10x Genomics) to produce single-cell gel bead-in-emulsions (GEMs). After reverse transcription, the GEMs were harvested, and the cDNAs were amplified and cleaned with the SPRIselect Reagent Kit (Beckman Coulter). Indexed sequencing libraries were constructed using the Chromium Single-Cell 3′ Library Kit (10x Genomics), for enzymatic fragmentation, end repair, A-tailing, adaptor ligation, ligation cleanup, sample index PCR, and PCR cleanup. The barcoded sequencing libraries were quantified by quantitative PCR, using the KAPA Library Quantification Kit (KAPA Biosystems). Sequencing libraries were loaded onto a NextSeq500 (Illumina), with a custom sequencing setting (26 bp for read 1 and 98 bp for read 2), to obtain a sequencing depth of approximately 80,000 reads per cell. Sequencing data were aligned to the human reference genome (GRCh38) and processed using the CellRanger 2.1.0 pipeline (10x Genomics).

Single-cell RNA-seq data preprocessing and clustering of chondrocyte cells

Raw gene expression matrices of all sample (including Blank sample, IL24 sample and IL48 sample) were combined in R (version 4.0.2) and converted to a Seurat object using the Seurat R package (version 3.2.0) [40]. Subsequently, all cells were removed that had either fewer than 500 UMIs, over 25000 or below 200 expressed genes, over 5% UMIs derived from mitochondrial genome, or log10 UMIs of per Gene lower than 0.8. Then we first used the NormalizeData function to normalize library size of each cell with default parameters and the FindVariableFeatures function to select the 3000 most variably expressed genes of each sample [41]. Next, we used the functions FindIntegration Anchors and IntegrateData as recommended by Seurat guidelines to integrated cells from three samples. This procedure created an integrated normalized dataset. We then scaled the integrated dataset using the ScaleData function. PCA analysis (Principal components analysis) was performed and the first 10 PCs, the number of principal components estimated by an Elbow plot, were used further to generated t-SNE dimensionality reductions of RunUMAP function. The clustering process was also performed based on these 10 PCs while graph-based clustering was run using FindNeighbors and FindClusters with a resolution of 0.5. [40, 42]. The corresponding marker genes were defined by FindAllMarkers function of Seurat package.

The correlation between clusters-special genes and functional gene sets and previously-identified marker genes of subpopulation of chondrocytes

Firstly, we select known gene sets that previous reported, for example, circadian rhythm (CRY1, CRY2, CLOCK, PER1, PER2, PER3, NR1D1, and BMAL1) [8], Microvasculature (ARAP3, ADCY4, ESAM, ERG, SLC43A3, SOX7, PTPRB, etc.) [43], ECM (MMP9, COL3A1, COL1A2, COL1A1, SOX5, PRG4, etc.), Proliferation (BAG1, ESR1, FOXA1, GPR160, NAT1, MAPT, MLPH, and PGR), and Inflammation (IFNG, IL5, IL6, JUN, NFKB1, STAT1, TGFB1, TNF, etc.) [44]. Detailed functional gene programs are listed in Supplementary Table 1. We calculated the average expression of each gene set in the total expression matrix. Then, the expression of corresponding top10 marker genes of each cluster defined by FindAllMarkers function of Seurat packages were also calculated. Correlations between the expression level of known gene sets and the expression level of the top 10 cluster-specific expressed genes were estimated by Pearson’s correlation coefficient. The markers of previously-identified chondrocyte subpopulations offered by previous study were also used to measure the similarity with the cell clusters by the same method The relevant code is shown in Supplementary File 1.

Single-cell trajectory reconstruction and analysis

MONOCLE (version 2.6.4) was applied to construct single-cell pseudo-time trajectories [45]. Briefly, a set of ordering genes, identified as being differentially expressed between clusters by Seurat, were selected and then passed to Monocle’s reversed graph embedding tool, which performed machine learning to learn a parsimonious principal graph. Dimensional reduction and cell ordering were performed using the DDRTree method and orderCells function with default parameters. The integrated high-dimensional expression profiles were reduced to a low-dimensional space, onto which each cell was projected and ordered in a trajectory. As branch points on the trajectory, the cells were divided into different “states,” and branch-dependent genes were then investigated by GO analysis using Metascape (https://metascape.org/) and Go database and Gene set variation analysis in the GSVA package (version 1.36.2).

Data available

The row data of single-cell RNA-seq are available at the NCBI’s Gene Expression Omnibus (GEO) data repository with the SRA accession ID PRJNA642343.