Article Text
PDF
Abstract
The cellular composition and underlying spatiotemporal transformation processes of subchondral bone in osteoarthritis (OA) remain unknown. Herein, various cell subsets from tibial plateau of patients with OA are identified, and the mechanism of subchondral microstructure alteration is elaborated using single-cell RNA sequencing technique. We identified two novel endothelial cell (EC) populations characterised by either exosome synthesis and inflammation response or vascular function and angiogenesis. Three osteoblast (OB) subtypes are introduced, separately related to vascularisation, matrix manufacturing and matrix mineralisation. The distinct roles and functions of these novel phenotypes in OA development are further discussed as well as interaction network between these subpopulations. The variation tendency of each population is testified in a destabilisation of the medial meniscus mouse model. The identification of cell types demonstrates a novel taxonomy and mechanism for ECs and OBs inside subchondral bone area provides new insights into the physiological and pathological behaviours of subchondral bone in OA pathogenesis.
- Osteoarthritis
- Therapeutics
- Biological Therapy
Data availability statement
Data are available in a public, open access repository. The single-cell RNA-seq data, quality control information and cluster information are available at the NCBI’s Gene Expression Omnibus (GEO) data repository with the accession ID GSE196678.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See: http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ON THIS TOPIC
Subchondral bone disorders during OA progress have been well recognized by advanced reports, but the specific cellular constituents and crosstalk are still unknown.
WHAT THIS STUDY ADDS
This study introduces cellular constituents and crosstalk inside subchondral bone from OA patients using single-cell sequencing technology, and verifies the findings in DMM mice model.
KDR- endothelial cells characterized by inflammatory response and KDR+ endothelial cells related with angiogenesis are identified, indicating aberrant vascularization inside OA subchondral bone.
MGST1+ EnOBs, GFPT2+ StOBs and WIF1+ MinOBs related with vascularization, matrix manufacturing and mineralization are pictured, indicating complicated functional changes of osteoblasts during OA progress.
HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE AND/OR POLICY
We raise new proof of subchondral bone based therapies for clinical research, and novel cellular/molecular targets for OA intervention.
Introduction
Osteoarthritis (OA) is an insidiously progressive, high-cost and poorly prognostic joint disease, affecting approximately 250 million patients worldwide.1 The most common clinical manifestations of OA are chronic cartilage degeneration, subchondral bone microstructure alteration, osteophyte formation and intractable joint pain.2–4 Fundamental research targeting articular cartilage destruction or cell senescence revealed considerable therapeutic potential of cartilage repair methods; however, the clinical trials have failed to varying degrees in recent decades.5–9 Accumulating evidence suggests that pathological alterations inside the subchondral bone are responsible for chondrocyte reduction and matrix degradation.10 11 The turnover rate of subchondral bone remains relatively low under normal circumstances and is accelerated by multiple factors in OA status, including mechanical and biological factors. The uncontrolled bone remodelling in the subchondral bone results in subsequent changes, including hypervascularisation, hyperpathia, abnormal mechanical support and cartilage destruction.12 Currently, both physiological and pathological behaviours of the subchondral bone have been valued in advanced research targeting OA therapies. However, the composition of subchondral bone cell types in patients with OA and the underlying spatiotemporal transformation processes remain unknown.
Mesenchymal stromal cells, osteoblasts (OBs), osteoclasts, endothelial cells (ECs) and immune cells are delicately orchestrated by various biological and mechanical factors in the local microenvironment of the subchondral bone.12 Under abnormal loading conditions, the activated form of transforming growth factor-β (TGF-β) is released from the bone matrix, resulting in aberrant vascularisation and osteogenesis.13 Hypertrophic chondrocytes also participate in this process as the major source of VEGF, coupled vessel invasion, cartilage remodelling and ossification.14 Interestingly, ECs recruited by multiple biological agents are the major driving forces of cartilage matrix degradation, bone elongation and remodelling.15 Moreover, OB-derived VEGF participates in the delicate bone–vessel crosstalk;16 however, the specific communication mode inside the OA subchondral bone remains unknown. OBs have multiple functions, including angiogenesis promotion, matrix manufacturing and mineralisation.17 Simultaneously, the participation of various immune cells leads to aggravated inflammation and subchondral bone disorders.18 This indicates the importance of subchondral bone cells in OA progression, and the multifaceted nature of ECs and OBs suggests that they are comprised of diverse subpopulations. Nevertheless, considering the high phenotypic heterogeneity and limited understanding of biomarkers, the isolation and definition of EC and OB subtypes in human subchondral bone remain unclear.
To reveal the cellular interactions involved in the OA subchondral environment, single-cell RNA sequencing (scRNA-seq) was performed on tibial subchondral bone samples from patients undergoing total knee arthroplasty. Here, the scRNA-seq technique was used to map a general census of subchondral bone cells from both normal and OA sites, determine the genetic characteristics of these cell subgroups and further analyse their potential differentiation relationships to characterise specific cell types. Finally, we investigated the cell–cell interaction network between EC and OB subpopulations. These results expand our understanding of the heterogeneity between patients and provide a theoretical basis for personalised OA therapies.
Methods
Human subchondral bone samples were collected during total knee arthroplasty operations. Bone samples were cut into 1–2 mm pieces and then digested in 0.2% type II collagenase (17101015, Thermo Fisher Scientific, Massachusetts) for 2 hours. Cells were then collected after red blood cell lysis.
Raw data processing and quality control
The Cell Ranger software pipeline (V.3.1.0) provided by 10× Genomics was used to demultiplex cellular barcodes, map reads to the genome and transcriptome using the STAR aligner, and down-sample reads as required to generate normalised aggregate data across samples, producing a matrix of gene counts versus cells. We processed the unique molecular identifier (UMI) count matrix using the R package Seurat (V.3.0). To remove low-quality cells and likely multiplet captures, which is a major concern in microdroplet-based experiments, we apply a criteria to filter out cells with UMI/gene numbers out of the limit of mean value ±2-fold of SD assuming a Guassian distribution of each cells’ UMI/gene numbers. Following visual inspection of the distribution of cells by the fraction of mitochondrial genes expressed, we further discarded low-quality cells where a certain percentage of counts belonged to mitochondrial genes. Library size normalisation was performed in Seurat on the filtered matrix to obtain the normalised count.
Top variable genes across single cells were identified using the method described in Macosko et al. Briefly, the average expression and dispersion were calculated for each gene, genes were subsequently placed into several bins based on expression. Principal component analysis was performed to reduce the dimensionality on the log-transformed gene-barcode matrices of top variable genes. Cells were clustered based on a graph-based clustering approach and were visualised in two dimension using t-stochastic neighbor embedding (tSNE). Likelihood ratio test that simultaneously test for changes in mean expression and in the percentage of expressed cells was used to identify significantly differently expressed genes between clusters. Here, we use the R package SingleR, a novel computational method for unbiased cell-type recognition of scRNA-seq to infer the cell of origin of each of the single cells independently and identify cell types.
Differentially expressed genes (DEGs) were identified using the Seurat package. P value <0.05 and |log2foldchange|>1 (or |log2foldchange|>0.58) were set as the threshold for significantly differential expression. Gene Ontology (GO) enrichment and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs were, respectively, performed using R based on the hypergeometric distribution.
Pseudotime analysis
Pseudotime analysis was performed with Monocle2 to determine the dramatic translational relationships among cell types and clusters. Further detection with the Monocle2 plot_pseudotime_heatmap function revealed the key role of a series of genes in the differentiation progress. Significantly changed genes were identified by the differential GeneTest function in Monocle2 with a q value<0.01.
Cell–cell communication analysis with CellPhoneDB 2
CellPhoneDB 2 is a Python-based computational analysis tool developed by Roser Vento-Tormo et al, which enables analysis of cell–cell communication at the molecular level. A website version was also provided for analysis of a relatively small data set (http://www.cellphonedb.org/). As described above, 606 single cells that were clustered into 5 cell types were investigated using the software to determine interaction networks. Interaction pairs including chemokines, ephrin receptor family, NOTCH family, cytokines and integrins and have p values <0.05 returned by CellPhoneDB were selected for the evaluation of relationships between cell types.
SCENIC analysis
Single-Cell rEgulatory Network Inference and Clustering (SCENIC) is a new computational method used in the construction of regulatory networks and in the identification of different cell states from scRNA-seq data. To measure the difference between cell clusters based on transcription factors (TFs) or their target genes, SCENIC was performed on all single cells, and the preferentially expressed regulons were calculated by the Limma package. Only regulons significantly upregulated or downregulated in at least one cluster, with adjusted p value <0.05, were involved in further analysis.
Animal models
Destabilisation of the medial meniscus (DMM) or sham operations were conducted bilaterally on 6-week-old male C57B6J mice (n=5 in each group). Samples were harvested at 0, 2, 4 and 8 weeks after surgery. Mice of 0 weeks were conducted with sham operation. Hearts were perfused with phosphate buffer solution (PBS) and 4% paraformaldehyde (PFA) successively in order to fix the antigens. Bilateral knee joints were harvested and fixed in 4% PFA for 24 hours, then went through micro-CT scan (60kV, 50μA, 10 um pixel).
Pathological staining
Samples were decalcified in 10% EDTA for 2 weeks, paraffin sections of 6 µm were then prepared for subsequent experiments: hematoxylin-eosin (HE), safranin O and fast green (SOFG) and immunofluorescence (IF) staining. HE and SOFG staining experiments were conducted under manufacturer’s instructions (Beyotime, C0105M; Solarbio, G1371). Antibodies used during IF staining: Rat anti Mouse PECAM1 (Thermo fisher, 1 40 311–82); Rabbit anti Mouse KDR (Abcam, ab11805); Rat anti Mouse OCN (Takara, M188); Rabbit anti Mouse MGST1 (Abcam, ab131059); Rabbit anti Mouse GFPT2 (Abcam, ab190966); Rabbit anti Mouse WIF1 (Santa, sc-373780); Goat anti Rat secondary antibody (Abcam, ba150165); Goat anti Rabbit secondary antibody (Abcam, ab150080). Nuclei were marked by DAPI (Beyotime, C1002).
Statistical analysis
Statistical calculations were performed using R package and GraphPad Prism (V.9.3). The results of microCT and staining are presented using line charts. The two-way Analysis of Variance (ANOVA) was applied to identify differences between groups in statistical graphs. Results were considered statistically significant when p<0.05.
Results
Single-cell profiling of human OA subchondral bone cells
To identify the cellular constitution of subchondral bone cells in human OA, we isolated human OA subchondral bone cells obtained from both lateral and medial tibial plateaus of two patients undergoing knee arthroplasty and profiled subchondral bone cells from different locations and patients (n=4) using scRNA-seq (figure 1 and online supplemental table S1). Through unbiased clustering of human subchondral bone cells, we found 10 clusters from patients with OA were identified, including T (11151), B (990), NK (5155), NKT (7538) and dendritic cells (DCs; 124), monocytes and macrophages (557), bone-related cells (864): ECs (246), mesenchymal stem cells (MSCs; 104) and OBs (360). According to animal experiments and clinical experience, OA involvement is more common and occurs earlier in the medial side of the tibial plateau. In this study, we selected patients with severe medial destruction and an almost healthy lateral plateau. These two patients were fully informed of their condition and chose total knee arthroplasty. Subchondral bone cells were divided into the control group (Ctrl; 13024) from the lateral tibial plateau and the OA group (13355) from the medial side (figure 1B). In total, 26 379 cells were retained for subsequent analysis after rigorous filtration (figure 1C,D; online supplemental figure S1A, S2A and table S2).
Supplemental material
Single-cell profiling of human osteoarthritis (OA) subchondral bone cells. (A) Schematic workflow of the experimental strategy. (B) X-ray photograph of a patient with knee OA (left) and corresponding cross-sectional anatomy of the subchondral bone and tibial plateau (right). (C,D) The t-stochastic neighbor embedding (tSNE) plots (left panel) and the sample origin (right panel) of 26 379 subchondral bone cells and 864 bone-associated cells. (E) Dot plots showing the distribution of each cell type in the control (Ctrl) and OA groups. Heatmap showing the pairwise correlations. (F) Cluster averaged log-normalised expression of the top 10 marker genes between the 10 cell types with stromal-related genes of interest annotated. Expression values are scaled per cluster. (G) Dot plot showing the expression of specific signatures in identified cell types in (F). The dot colour and size represent the mean expression and proportion of each cell population expressing genes, respectively. (H) Feature plots showing the expression of indicated markers for each cell type on the t-SNE map. EC, endothelial cell; MSC, mesenchymal stem cell; OB, osteoblast.
Next, the cell type distribution in the Ctrl and OA groups was analysed. In addition to the relationship between immune and myeloid cells, paired correlation analysis showed tight connections between ECs and MSCs, ECs and OBs and MSCs and OBs (figure 1E). The first 10 upregulated genes from these 10 clusters were used to create a heatmap (figure 1F). Representative markers for T, B, NK and NKT cells, DCs, monocytes, macrophages, ECs, MSCs and OBs were revealed (figure 1G). Specifically, the following clusters were identified: (1) B cells (expressing CD79A, BANK1 and MS4A1), (2) NK cells (expressing GZMB and NKG7), (3) NKT cells (expressing NKG7 and CD3D), (4) T cells (expressing CD3D and CD3G), (5) DCs (expressing LILRA4 and PTCRA), (6) monocytes (expressing CSTA and FCN1), (7) macrophages (expressing CD14, CD68, CSF1R, C1QC and F13A1), (8) ECs (expressing PECAM1 and CLDN5), (9) MSCs (expressing MCAM) and (10) OBs (expressing RUNX2 and cadherin 11 (CDH11); figure 1 and online supplemental figures S1B and S2B–D).
Identification of bone-related cell populations in human OA subchondral bone
Abnormal angiogenesis, subchondral bone remodelling and sensory innervation are well recognised during early stage of OA and might cause cartilage destruction and pain directly or indirectly.13 Therefore, OBs, osteoclasts, ECs and neuronal cells were the focus of our research, rather than immune cells. To define the cell subpopulation and identify genome-wide gene expression patterns, bone-related cells were clustered to produce eight clusters (figure 2). Next, to explore the potential transformation between different cell types and visually depict the differentiation paths, the Monocle method was used to determine the pseudotemporal order between cell types (figure 2B,C). Our analysis clearly identified and verified three major groups of differentiated cell types: ECs (PECAM1+), MSCs (MCAM+) and OBs (RUNX2+/CDH11+; figure 2D).
Identification of bone related cell populations in human osteoarthritis (OA) subchondral bone. (A) t-SNE plots of bone associated cells coloured by cluster. (B–D) Pseudotime trajectory plot showing differentiated cell types (endothelial cells (ECs), mesenchymal stem cells, and osteoblasts (OBs) at the end of the branches. Dots along the trajectory lines represent the status of the cells transitioning toward differentiated cell types. (E) Heatmap revealing the scaled expression of differentially expressed genes (DEGs) for each cluster defined in (A). (F–H) Violin plots showing expression levels of indicated markers for eight clusters. (I–J) Single-cell regulatory network inference and clustering analysis showing distinct regulons in eight clusters. The heatmap shows only the regulons with significant differences. MSC, mesenchymal stem cell.
Among the eight clusters, cluster1 expressed markers of multiple cell types, such as CTSK, RGS10 and SPP1, suggesting that cluster1 contains osteoclasts, nerve cells and OBs (online supplemental figure 3). Cluster1 is a heterogeneous cell cluster; therefore, we only analysed OBs, ECs and MSCs. To determine the characteristics of each cell cluster by analysing differential gene transcript expression patterns, a differentially expressed feature analysis was performed using the scRNA-seq data set and all cell clusters were compared with one another. We discovered 78 DEGs that best divided bone-related cells in subchondral bone into eight subclusters (figure 2E). Next, the original sample information and expression levels of the indicated markers were combined to determine the cell identity of each cluster, and their biological functions were analysed via regulons CSI correlation heatmap of the coexpression between TFs and potential target genes. Seven major cell clusters were identified: precursor ECs, pre-ECs (C2CD4B+/B3GNT5+); ECs (VWF+/KDR+); endothelial OBs (EnOBs, ABCA10+/microsomal glutathione S-transferase 1 (MGST1)+); stromal OBs (StOBs, PTGS2+/glutamine-fructose-6-phosphate transaminase 2 (GFPT2)+); mineralised OBs (MinOBs, WNT inhibitory factor 1 (WIF1)+/NDNF+) and two MSC subpopulations (figure 2F–H). Through SCENIC analysis, we discovered that pre-ECs and ECs exhibit activated proangiogenesis regulons, such as SMAD Family Member 1 (SMAD1), ETS transcription facto (ERG) and ETS1, and that the regulators have higher activity in ECs than in pre-ECs. Furthermore, OB subpopulations exhibited similar activated TFs; however, the regulatory activities of TFs are different (figure 2I,J).
Identification of pre-ECs and ECs
In addition to the well-known differences between arteries, capillaries and veins, ECs are highly heterogeneous and acquire specialised functional properties in the local microenvironment. The articular cartilage is constantly maintained in a low-oxygen environment. However, due to the high metabolic requirements of OBs, blood vessels are required to provide sufficient oxygen. The cells are relatively hypoxic during osteogenesis. Osteogenesis and nearby ECs may increase HIF-1α expression. Upregulation of HIF-1α activity in hypoxic tissues leads to increased VEGF expression and promotes angiogenesis. Using IF and flow cytometry to identify the EC subpopulation in bone, Kusumbe et al proposed the following terminology for bone microvessels: H-type for the small PECAM1hi/Emcnhi subset and L-type for the PECAM1lo/Emcnlo sinusoidal vessels.19 As previously mentioned, EC subpopulations were divided into pre-ECs and ECs, and we discovered certain differences between them (figure 2E,F). To investigate the distinct features of pre-ECs and ECs, we identified DEGs between them (figure 3A). The GO and KEGG were analysed with these DEGs to show pre-EC and EC characteristics. Notably, pre-ECs were enriched for extracellular exosomes, interleukin-mediated signalling pathways and ribosomes, whereas ECs were enriched for vasculogenesis, angiogenesis, EC migration and signalling pathway regulation, such as VEGF, Rap1, PI3K/Akt, Ras and MAPK signalling pathways (online supplemental figure S4A–D). Then we compared characteristic genes of these two clusters Pre-EC identification markers have diverse functions, including genes related to ribosome synthesis, exosome synthesis and inflammation, including RPL17, HNRNPF, RABA5 and C2CD4B, whereas ECs are primarily enriched in genes related to angiogenesis, such as VWF, KDR, TIE1 and CDH5 (figure 3B). Moreover, angiogenesis-related EMCN, PECAM1, EGFL7, ENG and KDR were all upregulated during the differentiation of pre-ECs to ECs, while exocrine-related DDIT3 and RAB5A73 and inflammation-related CCL2 and C2CD4B74 were downregulated (figure 3C).
Identification of precursor ECs (pre-ECs) and ECs. (A) Heatmap of DEGs between different ECs. (B) Violin plots showing the expression levels of the specific representative genes marking pre-ECs and ECs. (C) Heatmap showing upregulation or downregulation of vascular markers, exosomes, and ribosomal markers in the differentiation process. (D) tSNE plots of the expression levels of transcription factors (TFs; up) and area under the curve scores (down). (E) Kyoto Encyclopaedia of Genes and Genomes pathway enrichment between the OA and Ctrl groups in pre-ECs and ECs. (F) Gene Set Enrichment Analysis (GSEA) pathway enrichment between the OA and Ctrl groups in pre-ECs and ECs. EC, endothelial cell; DEGs, differentially expressed genes; OA, osteoarthritis.
Next, the essential motifs of the two EC subpopulations were identified using SCENIC analysis. As the specific motifs of pre-ECs, ATF3 and CCAAT Enhancer Binding Protein Delta (CEBPD) are essential in the transcriptional regulation of inflammation, whereas ERG and SMAD1 motifs, which are closely related to angiogenesis, are highly activated in ECs. SOX17 and JunD are TFs that are widely expressed in pre-ECs and ECs (figure 3D and online supplemental figure S5A,B). In a previous study, SMAD1 sprouted angiogenesis in human embryonic stem cell-derived ECs.20 ERG is an essential regulator of angiogenesis and vascular stability through Wnt signalling.21 These results may help us to identify pre-ECs and ECs and expand our understanding of the novel function of subchondral EC subsets in OA.
During OA progression, the EC cluster was a substantially increased cell population (figure 1E). To investigate the distinct features of the OA and Ctrl groups, we identified DEGs between these two groups in pre-ECs and ECs (online supplemental figure S6A,B). To analyse the identified features of pre-EC and EC clusters from the OA group, we analysed the differences using GO, KEGG and Gene Set Enrichment analyses (GSEA). Notably, compared with that of the Ctrl group, pre-ECs of the OA group showed stronger protein synthesis, inflammation-related pathways and responses, whereas ECs were enriched for angiogenesis-promoting functions, such as blood vessel development, EC differentiation and platelet-derived growth factor binding (figure 3E,F and online supplemental figure S6C,D).
Determining the relationships among EnOBs, StOBs, and MinOBs
According to Rutkovskiy et al, OBs undergo a three-stage differentiation: stage 1, the cells continue to proliferate; stage 2, they start differentiating, while maturing the extracellular matrix (ECM) with alkaline phosphatase (ALP) and collagen and stage 3, the matrix mineralisation and mineral deposits increase.22 According to the SCENIC analysis (figure 2I), we found that the TF types and transcriptional activities were different among EnOBs, StOBs and MinOBs, which means that these clusters may have different cellular functions. By identifying DEGs between EnOBs, StOBs and MinOBs, the differences between OB subpopulations were further verified (figure 4A). To explain the specific characteristics of these three OB populations, we analysed the differences through GO and KEGG analyses. EnOBs are related to EC migration, VEGF binding and the PDGFR-β signalling pathway and express NRP1, PDGFRB and VCAM, suggesting that this cluster may have potentially affect angiogenesis. StOBs were enriched for collagen and fibre-related biological processes, such as collagen fibril organisation, fibronectin binding, and ECM binding. MinOBs distinctively expressed an ossification and bone mineralisation biological process gene signature (figure 4B and online supplemental figure S7A–C). Next, representative candidate markers among EnOBs, StOBs and MinOBs (figure 4C) were identified, including TFs (figure 4F).
Determining the relationships among endothelial osteoblasts (EnOBs), stromal OBs (StOBs) and mineralised OBs (MinOBs). (A) Heatmap showing Z score scaled expression levels of DEGs for EnOB, StOB, and MinOB populations. (B) Heatmap showing the differences in enriched gene ontology (GO) functions of upregulated genes in different OB subsets. (C) Boxplots showing the expression levels of representative candidate marker genes specifically expressed in different subsets. (D,E) Monocle pseudospace trajectory revealing the progression of OB lineage in subchondral bone coloured according to cluster. Monocle pseudotime trajectory revealing the progression of EnOBs, StOBs and MinOBs. (F) Pseudotemporal expression dynamics of TFs in EnOBs, StOBs, and MinOBs. All single cells in the EnOB, StOB and MinOB cell lineage are ordered based on pseudotime. (G) GO functions enrichment analysis of OA versus Ctrl upregulated genes in EnOBs, StOBs and MinOBs. (H) GSEA showing enrichment of pathways between the OA and Ctrl groups in EnOBs, StOBs and MinOBs. ECM, extracellular matrix; DEGs, differentially expressed genes; GSEA, Gene Set Enrichment Analysis; OA, osteoarthritis; TFs, transcription factors.
Next, the Monocle method was applied to depict the pseudotemporal sequence of potential differentiation pathways among cell types. According to the pseudotime trajectory axis, we suggest that StOB is an intermediate state representing the state between EnOBs and MinOBs (figure 4D,E). The pseudotemporal expression dynamics of representative candidate markers and TFs also marked the progression from EnOBs to StOBs to MinOBs. Through SCENIC analysis, we found that TFs upregulated by EnOBs were related to angiogenesis, such as MECOM, ERG and XBP1, and TFs related to collagen and fibre in StOBs, including FOSL2, ERG1 and ATF4, and VDR and FOXC2 in MinOBs are related to mineralisation (figure 4F and online supplemental figure S8A–C). Taken together, these data reveal the relationships and potential functions of EnOBs, StOBs and MinOBs.
The number of OBs increased in abundance during OA progression, more so than in ECs (figure 1E). To investigate the distinct features of normal and diseased cells inside the subchondral zone, we identified DEGs between the OA and Ctrl groups of these three OB clusters, respectively (online supplemental figure S9A–C). We then analysed the differences by GSEA and GO analysis. Notably, compared with that of the Ctrl group, the OA group of EnOBs showed stronger angiogenesis and wound healing biological processes; the OA group of StOBs was enriched for ECM binding and collagen fibril organisation, and the OA group of MinOBs was enriched for response of metal ions such as cadmium, copper and zinc (figure 4G,H). Regarding these ionic reactions, copper ions in biological materials promote bone formation,23 zinc increases OB activity and collagen synthesis24 and cadmium promotes OB differentiation.25 These results indicate that OBs in the OA group have stronger osteogenic effects and we will further verify the functional heterogeneity of the three OB subpopulations.
Vascular EC and OB subpopulation interaction
Through the interaction of H-type blood vessels and various cytokines in bone metabolism, angiogenesis and bone formation are precisely coupled.19 MSCs are chemically recruited by ECs to promote osteogenesis.26 Furthermore, the crucial signalling pathways in MSCs coupled with ECs include the TGF-β, PDGF-PDGFR, angiopoietin, Notch and FAK signalling pathways. To further explore the key signalling pathway that couples OBs and ECs, we studied the cell–cell interaction network between the identified clusters. Considering the outcomes of GSEA and GO enrichment analyses and the characteristics of the genes and TFs, we chose to analyse the interaction pairs, including chemokines, ephrin receptor family, NOTCH family, cytokines and integrins. We found that ECs were the predominant cell population interacting with the OB subpopulation, and pairwise correlation analysis revealed that ECs are more closely related than pre-ECs to OB subpopulations (figure 5A and B).
Vascular EC and OB subtype interaction. (A) Pearson correlation analysis of two clusters of ECs subsets and three clusters of OBs. (B) CellPhoneDB analysis showing the number of ligand-receptor interactions between EC and OB subpopulations. Circos plots showing ligand-receptor pairs of cytokines (C), growth factors (D), and integrin (E) between EC subpopulations. Bubble plots showing ligand-receptor pairs of cytokines (F), growth factors (G), and integrin (H) between OB subpopulations. EC, endothelial cell; EnOBs, endothelial OBs; MinOBs, mineralised OBs; OB, osteoblast; StOBs, stromal OBs.
Compared with pre-ECs, ECs express a larger number of membrane receptors, fibronectin and collagen and secrete a larger number of angiocrine factors. EC analysis revealed that ECs exhibited abundant expression of multiple membrane receptors for ligands important for vascular development, including NOTCH1, NOTCH4, VEGF receptors (KDR, FLT1, FLT4, NRP1, NRP2), TGFβ receptors (TGFBR3), ephrin B receptor (EPHB4) and tyrosine kinase receptor (TEK), which bind to JAG1, the VEGF family, PGF, ANGPT1 and TGFβ 1 ligand secreted by OB to promote angiogenesis (figure 5C–E and online supplemental figure S10A–C). These results indicate that ECs are a mature EC subgroup with angiogenic function at the transcriptional level and are mainly coupled with OBs. Notably, pre-ECs hardly secrete any ligands, as the results showed, however, considering the enrichment analysis data of pre-ECs (online supplemental figure S4A), we hypothesise that the function of pre-ECs is achieved by secreting exosomes, which requires further investigation.
There is little difference between chemokine interaction pairs in OBs, and they all express CXCR4 to promote OB proliferation and differentiation; however, only MinOBs do not secrete CXCL12.27 At the end of osteogenic differentiation, CXCL12 is downregulated.28 We found that the overall expression of the ephrin receptor and NOTCH family members decreased gradually from EnOBs to StOBs to MinOBs. The ephrin receptor29 and NOTCH families30 promote OB proliferation, and the lack of the bone system JAG1 leads to mature OB proliferation, which is manifested as an increase in the rate of mineral deposition (figure 5F). Additionally, we found that these three OB subgroups are affected by bone morphogenetic proteins, which have a strong positive effect on bone formation.31 PGF, PDGF and VEGF families secreted by ECs and OBs could interact with the receptors highly expressed on the surface of EnOBs, including NRP1, PDGFRA and PDGFRB. These cytokines are related to angiogenesis, and PDGF induces OB proliferation via the ERK signalling pathway.32 We also found that TGFβ1 was produced by all three types of OBs, but only bound with TGFβ receptors on EnOBs to play a role in promoting OB proliferation,33 inducing VEGF secretion34 and inhibiting mineralisation function.35 As a non-collagenous protein in the bone ECM that is recognised to regulate bone formation and mineralisation, osteopontin (OPN/SPP1) is expressed and released in the integrin interaction pair by MinOBs only (figure 5G,H). Compared with that of EnOBs and MinOBs, StOBs expressed the highest fibrin and collagen levels, while MinOBs expressed the lowest (figure 5H). In summary, the above data further validated our understanding of the role of these five major cell clusters in bone-associated cells and the interaction between ECs and OBs.
Pathological identification of subpopulations
In order to further verify our sequence results, we conducted DMM surgery on 6-week-old C57B6J mice, simulating patient conditions before arthroplasty surgery (figure 6A). The tibial subchondral bone volume in OA mice showed significant changes after surgery, as shown by microCT, safranin O, fast green staining and H&E staining (figure 6B). The total tissue volume of subchondral bone decreased at 2 weeks and increased at 4 weeks postsurgery, and the subchondral bone structure densified at 8 weeks (figure 6C). Microstructure disruption was indicated by aberrant subchondral bone plate thickness and trabecular pattern factor, according to the microCT calculation (figure 6C). Osteoarthritis Research Society International scores indicated that cartilage degeneration was significant at 4 weeks and deteriorated at 8 weeks postsurgery (figure 6D).
Pathological identification of DMM mice model. (A) Schematic illustration of the experimental process. (B) Top row, representative three-dimension image of subchondral bone from tibial medial plateau at 0, 2, 4 and 8 weeks after sham or DMM surgery. Middle rows, safranin O and fast green stain of tibial medial plateau at the same checkpoint. Bottom rows, H&E staining of knee joint, cartilage and subchondral bone. (C) Quantitative analysis of bone volume tissue/total tissue volume (BV/TV), subchondral bone plate thickness (SBP. Th) and trabecular bone pattern factor (Tb.Pf) in medial plateau calculated from micro-CT results. (D) Osteoarthritis Research Society International scores after DMM operation are shown bottom right. n=10 in each group. Scale bar, 100 µm. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001. Statistical significance was shown by two-way ANOVA. Data were presented as the mean±SD. DMM, destabilisation of the medial meniscus; HC, hyaline cartilage; CC, calcified cartilage; SBP, subchondral bone plate.
An IF test was performed to characterise cell development and trajectory during OA pathology. The ratio of ECs among the total PECAM1-positive cells increased continuously after surgery (figure 7A). Specifically, PECAM1-positive ECs were rather rare and were near the trabecular before intervention, indicating a relatively low angiogenesis rate. An important proportion of PECAM1-positive cells was KDR-negative or pre-ECs, in the sham group (66.4%±9.44%) and 0w samples in the DMM group (60.69%±8.77%). They are characterised by genes coded for ribosome synthesis, extracellular vesicles synthesis and inflammation, indicating hypermetabolism and proinflammatory status inside the subchondral bone. ECs grew more after 4 weeks and the ratio of KDR-positive, or ECs with relatively high angiogenesis trends, increased significantly and reached 82.97%±8.01%. With severe subchondral sclerosis progression, the majority (88.90%±4.62%) of ECs became pre-ECs, and the total number of pre-ECs and ECs decreased dramatically in the limited space (figure 7B).
Pathological identification of EC and OB subpopulations. (A) Immunofluorescence staining of PECAM1, KDR, white arrow shows co-positive area of corresponding markers. (B) Percentage of PECAM1 +KDR+ ECs in total PECAM1 +cells. (C) Immunofluorescence staining of OCN, MGST1, GFPT2 and WIF1, white arrow shows co-positive area of corresponding markers. (D) Percentage of MGST1 +EnOBs, GFPT2 +StOBs and WIF1 +MinOBs in OCN+cells, n=9 in each group. Scale bar, 100 µm. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001. Statistical significance was shown by two-way ANOVA. Data were presented as the mean±SD. ANOVA, analysis of variance; EC, endothelial cell; EnOBs, endothelial OBs; MinOBs, mineralised OBs; OB, osteoblast; StOBs, stromal OBs; SHAM, sham operation; DMM, destabilization of the medial meniscus.
Next, we analysed the OB subpopulation marked by osteocalcin (OCN), MGST1, GFPT2 and WIF1. Consistent with the previous results, the total number of OCN-positive cells increased during the first 4 weeks and decreased at 8 weeks postsurgery (figure 7C). Analogously, EnOBs characterised by angiogenesis-related genes increased dramatically 2 weeks post-DMM surgery and decreased at 8 weeks (figure 7D, top). StOBs capable of ECM binding and collagen fibril organisation increased continuously during the first 4 weeks post-trauma and decreased at 8 weeks (figure 7D, middle), probably because of the calcification requirement. MinOBs, which were closely related to metal ion and biomineralisation, increased dramatically at 8 weeks postsurgery and accounted for approximately 48% of the total OCN-positive cells (figure 7D, bottom). Taken together, these data indicate chronological changes in OB subgroups and mapping the over-time changing bio-function of OBs in subchondral bone during post-traumatic OA.
Discussion
OA is one of the most common degenerative diseases that cause disability in older adults. An epidemiological study by Tang et al showed that 8.1% of the adult population had clinically significant OA of the knee or hip.36 Moreover, OA consumes a substantial amount of healthcare resources, primarily owing to the joint replacement surgery costs for advanced OA. Increasing economic pressure, an ageing society and the obesity epidemic emphasise the need for new strategies for the diagnosis and intervention at early-stage OA.37–39 Increasing evidence suggests that the appearance of subchondral bone lesions occurs earlier than cartilage degeneration, and the pathological alterations of subchondral bone play an important role in OA development. During the initial phase of OA, the bone turnover rate beneath the articular cartilage was upregulated, and vascular invasion took place bottom-up from the subchondral, the tidemark and ultimately into the cartilage. Since the exact role of subchondral bone during OA initiation and progression remains unclear, and the specific cell markers are lacking, it is urgent to unveil the internal state of subchondral bone cells in OA pathogenesis. Here, we used comprehensive gene expression profiling to reveal the cell types that make up the subchondral bone microenvironment at single-cell resolution and also novel cell markers and characteristics to verify each hypothetical subchondral bone cell cluster.
We identified 10 different cell types in the human OA subchondral bone microenvironment. Notably, there were more OBs and ECs in the OA group than in the Ctrl group. This finding validates the characteristics of increased bone formation and angiogenesis in OA subchondral bone.40 In addition to the empirically inferred subchondral bone cell types, we identified new subtypes of bone-associated cells and new markers of bone-associated cell populations based on scRNA-seq analysis. Based on the expression of TFs and markers of the new subtypes, we suggest that EC subtypes and OB subtypes perform different biological functions.
In the process of OA cartilage erosion, compared with the top-down vessel invasion originating from synovial tissue or synovium, bottom-up vascularisation from subchondral bone plays a larger role.41 42 In the process of bone elongation, ECs constantly erode the cartilage matrix, thus creating space for osteogenesis. Kusumbe et al identified two types of special blood vessel subtypes based on the expression strength of PECAM1 and EMCN, namely, H-type (PECAM1hiEMCNhi) and L-type blood vessels (PECAM1loEMCNlo).19 In the present study, we also identified two EC phenotypes: pre-ECs and ECs. We found that PECAM1 and EMCN expression in ECs was upregulated, suggesting that it may have functions similar to those of H-type blood vessels. By comparing gene expression, TF activity and enrichment analysis, we suggest that the ‘ECs’ in this research are a type of ECs that promote angiogenesis and can also be coupled with OBs. In this study, we identified a new subset named pre-ECs, characterised by interleukin-mediated inflammation pathways and exosomes. A key element leading to the advancement of OA is the production of high inflammatory cytokine levels, and the proinflammatory cytokine interleukin 1β (IL-1β) is expressed in large quantities during OA. Elevated IL-1β levels are associated with tissue damage and reflect the severity of inflammation. During subchondral bone reconstruction, IL-1β may also play a role in promoting cartilage calcification (ossification) and cartilage degeneration.43 Yang et al showed that exosomes derived from vascular ECs promoted the progression of OA by promoting chondrocyte apoptosis.44 Therefore, we speculate that pre-ECs in OA may allow for bone formation through exosomes digesting cartilage and promote subchondral bone remodelling through IL-mediated inflammation. The characterisation of these two novel subsets improves our understanding of the characteristics and functions of ECs in the OA subchondral bone.
Subchondral bone sclerosis is characterised by an increase in bone volume due to an enhanced bone turnover rate.40 In the OA subchondral bone, we found three OB phenotypes: EnOBs, StOBs and MinOBs. We believe that EnOBs are enriched in angiogenesis-related pathways, StOBs are characterised by collagen and fibrosis, and MinOBs specifically express mineralisation-related markers. Here, we found that the three OB phenotypes have a differentiation path, from EnOBs to StOBs to MinOBs. The first stage of OB differentiation is characterised by cell proliferation. EnOBs are tightly correlated with vascularisation, including EC migration, VEGF binding and platelet-derived growth factor receptor β (PDGFRβ) signalling pathways. Angiogenesis directly promotes bone formation, and CDH11 cooperates with PDGFRβ to promote cell proliferation. In the next stage, OBs begin to differentiate and express collagen and ALP. StOBs are OBs characterised by OB differentiation and collagen fibril organisation and express high levels of fibulin-1 (FBLN1), collagen type XI alpha 1 (COL11A1) and PLOD2. FBLN1 is an important ECM protein that stabilises collagen and other ECM proteins. COL11A1, as one of the three alpha chains of type XI collagen, is critical for collagen fibre assembly and bone development. OA-related fibrosis is associated with elevated PLOD2 expression.45 In the last stage, OBs express markers of more bone sialoprotein, OPN/SPP1 and OCN (BGLAP), thereby inducing matrix mineralisation. In this study, we found that MinOBs specifically express the above-mentioned mature OB phenotypic markers, representing the biological processes of ossification, bone mineralisation and biomineral tissue development. Additionally, MinOBs express the Wnt antagonist WIF1, which regulates the Wnt/β-catenin signalling pathway to reduce cell proliferation and promote mineralisation.46
The coupling of osteogenesis and angiogenesis causes increased bone mineral density and significant microstructural changes in the subchondral bone. Through the analysis of cell-to-cell communication between ECs and OBs, we found that ECs with angiogenic function and upregulated PECAM1 and EMCN expression are coupled with OBs through NOTCH, VEGF, TGFβ receptors, EPHB4 and TEK. Although we highlighted the role of ECs in cell-to-cell communication instead of pre-ECs, the contribution of pre-ECs in the advancement of OA could not be neglected. We found that specific receptors on OBs, including EPH, NOTCH, BMPR, NRP, PDGFR, TGFβR, fibronectin, collagen and SPP1 (OPN), responded to signals derived from ECs, thus further elucidating the interaction between the osteogenic subgroups and ECs. Furthermore, EnOBs are related to angiogenesis and promote OB proliferation. StOBs express higher levels of fibrin and collagen. Moreover, OPN, which is directly related to mineralisation, is expressed only by MinOBs.
In conclusion, our scRNA-seq analysis results provided a clearer and more consistent definition of the cellular components of human subchondral bone in OA. Specifically, two novel populations of ECs and three subpopulations of OBs, as well as the intercellular interaction network between these subpopulations, were identified. Our analysis provides new insights into the physiological and pathological behaviours of subchondral bone in OA pathogenesis, which may contribute to novel therapeutic strategies in the future.
Data availability statement
Data are available in a public, open access repository. The single-cell RNA-seq data, quality control information and cluster information are available at the NCBI’s Gene Expression Omnibus (GEO) data repository with the accession ID GSE196678.
Ethics statements
Patient consent for publication
Ethics approval
This study involves human participants and was approved by Ethical Committee of Shanghai University (ECSHU 2021-146). Participants gave informed consent to participate in the study before taking part.
Acknowledgments
We thank OE biotech company (Shanghai, China) for the support of bioinformatics analysis.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
YH, JC, HL and SW contributed equally.
Contributors YH: conceptualisation, methodology, investigation, writing. XC: methodology, visualisation. JC: validation, writing original draft. SW: formal analysis, software. QZ: methodology, sample preparation. HL: sample preparation. HZ: visualisation. JG: investigation, data curation. LC: investigation, project administration, data curation. KX: supervision, funding acquisition, writing. JS: supervision, funding acquisition and guarantor.
Funding This work is supported by National Key R&D Program of China (2018YFC2001500), National Natural Science Foundation of China (82172098) and Shanghai Municipal Health Commission (202040372).
Competing interests None declared.
Provenance and peer review Not commissioned; externally peer reviewed.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.