Article Text

CD_99 G1 neutrophils modulate osteogenic differentiation of mesenchymal stem cells in the pathological process of ankylosing spondylitis
  1. Xinzhe Feng,
  2. Chen Wang,
  3. Boyao Ji,
  4. Junjie Qiao,
  5. Yihong Xu,
  6. Shanbang Zhu,
  7. Zhou Ji,
  8. Bole Zhou,
  9. Wenwen Tong,
  10. Weidong Xu
  1. Department of Joint Bone Disease Surgery, Changhai Hospital, Naval Medical University, Shanghai, China
  1. Correspondence to Professor Weidong Xu, Department of Joint Bone Disease Surgery, Changhai Hospital, Shanghai, Shanghai 200433, China; xuwdshanghai{at}126.com; Dr Wenwen Tong, Department of Joint Bone Disease Surgery, Changhai Hospital, Shanghai, Shanghai 200433, China; tww944057350{at}163.com

Abstract

Objectives This study aimed to identify the types and heterogeneity of cells within the spinal enthesis and investigate the underlying mechanisms of osteogenesis.

Methods Single-cell RNA sequencing was used to identify cell populations and their gene signatures in the spinal enthesis of five patients with ankylosing spondylitis (AS) and three healthy individuals. The transcriptomes of 40 065 single cells were profiled and divided into 7 clusters: neutrophils, monocytic cells, granulomonocytic progenitor_erythroblasts, T cells, B cells, plasma cells and stromal cells. Real-time quantitative PCR, immunofluorescence, flow cytometry, osteogenesis induction, alizarin red staining, immunohistochemistry, short hairpin RNA and H&E staining were applied to validate the bioinformatics analysis.

Results Pseudo-time analysis showed two differentiation directions of stromal cells from the mesenchymal stem cell subpopulation MSC-C2 to two Cxcl12-abundant-reticular (CAR) cell subsets, Osteo-CAR and Adipo-CAR, within which three transcription factors, C-JUN, C-FOS and CAVIN1, were highly expressed in AS and regulated the osteogenesis of mesenchymal stem cells. A novel subcluster of early-stage neutrophils, CD99_G1, was elevated in AS. The proinflammatory characteristics of monocyte dendritic cell progenitor—recombinant adiponectin receptor 2 monocytic cells were explored. Interactions between Adipo-CAR cells, CD99_G1 neutrophils and other cell types were mapped by identifying ligand–receptor pairs, revealing the recruitment characteristics of CD99_G1 neutrophils by Adipo-CAR cells and the pathogenesis of osteogenesis induced in AS.

Conclusions Our results revealed the dynamics of cell subpopulations, gene expression and intercellular interactions during AS pathogenesis. These findings provide new insights into the cellular and molecular mechanisms of osteogenesis and will benefit the development of novel therapeutic strategies.

  • Ankylosing Spondylitis
  • Inflammation
  • Autoimmune Diseases

Data availability statement

Data are available upon reasonable request. Not applicable.

http://creativecommons.org/licenses/by-nc/4.0/

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

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.

WHAT IS ALREADY KNOWN ON THIS TOPIC

  • Little is known regarding the underlying mechanisms of osteogenesis in ankylosing spondylitis (AS), a chronic inflammatory disease characterised by bony overgrowths of the axial spine.

WHAT THIS STUDY ADDS

  • On revealing the cell profiles of spinal enthesis in patients with AS, we found a close interaction between early-stage CD99_G1 neutrophils and Adipo-Cxcl12-abundant-reticular cells in the pathogenesis of AS. CD99_G1 neutrophils also express various secreted proteins contributing to osteogenesis. Additionally, MDP-ADIPOR2 monocyte subclusters participate in multiple pathological processes associated with AS, many of which are related to inflammation.

HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY

  • Our work contributes to a better understanding of the pathogenesis of AS and highlights novel targets for AS prevention or treatment.

Introduction

Ankylosing spondylitis (AS) is a chronic autoimmune inflammatory disease of the axial skeleton that triggers ossification and causes spinal deformities.1 2 It is estimated that the ratio of male to female in AS was 2-3:1, and AS characteristically affects young people, with a peak age between 20 and 30 years.3 Gradual functional spinal impairment makes AS one of the most common causes of severe disability worldwide, affecting millions of people. Currently, because of various biologics and non-steroidal anti-inflammatory drugs, patients suffering inflammatory pain caused by AS can achieve short-term relief, though ankylosis is nearly irreversible. The prevailing view is that ossification of the spine in AS begins at the enthesis where the tendon and ligament attach to the bone near the spine.4 As an immune-mediated inflammatory disorder, the role of immune cells is crucial, and some functions of immune cells have been discovered in pathological processes in AS.5–7 Inductive factors of AS and their underlying mechanisms have not been elucidated.

Single-cell RNA sequencing (scRNA-seq), a high-quality alternative to traditional transcriptome screening, has enabled the simultaneous characterisation of gene expression levels in thousands of individual cells with high throughput and high resolution.8 This method has delivered revolutionary insights into cellular dynamics and heterogeneity and has the potential to investigate cell-to-cell variations, identify possible subpopulations and discover possible rare cell types.9 Although scRNA-seq is increasingly being applied, its application in AS has been limited to peripheral blood mononuclear cells and synovium of knee joint,10–13 leaving the transcriptomic heterogeneity within the pathological tissue of spine largely unknown.

In this study, scRNA-seq was employed at the spinal ligament attachment site, a pathologically affected lesion in patients with late-stage AS. We aimed to determine the transcriptomic profiles, heterogeneity and subpopulation composition of cells from the spinal enthesis of five patients with AS and three healthy controls. We also aimed to uncover the underlying mechanisms by which AS induces mesenchymal stem cell (MSC) osteogenesis (online supplemental material 1).

Supplemental material

Results

Single-cell transcriptome landscape in spinal ligament attachment site of AS

To understand the microenvironment of spinal enthesis in AS at single-cell resolution, we collected and profiled spinal tissue samples from eight patients (five patients with AS and three healthy controls) (figure 1A, online supplemental table S1). A total of 40 065 cells were recovered according to the quality control criteria, of which 23 683 and 16 382 cells were from the AS and control groups, respectively. The number of genes and unique molecular identifiers (UMIs) detected per cell and the percentage of mitochondrial genes were all within the normal range (online supplemental tables S2 and S3). Unbiased clustering of all included cells yielded 21 clusters according to principal component analysis as visualised by uniform manifold approximation and projection analyses (figure 1B, online supplemental figure S1A) covering seven different cell types identified based on marker genes. These were neutrophils (S100A8+ S100A9+), monocytes (IL-1β+ LYZ+), granulomonocytic progenitor (GMP)_erythroblasts (GATA2+ CD34+), T cells (CD3E+ CD3D+), B cells (CD19+), plasma cells (JCHAIN+) and stromal cells (SPARC+ DCN+) (figure 1C, online supplemental figure S1B, online supplemental table S4).

Supplemental material

Supplemental material

Supplemental material

Supplemental material

Supplemental material

Supplemental material

Figure 1

Overview of the single-cell landscape for the spinal ligament attachment site in ankylosing spondylitis (AS) and healthy control patients. (A) Processing pipeline of samples from patients with AS and healthy controls. The single-cell suspension was separated from spine entheseal tissues and subsequently loaded onto the BD scRNA-seq platform. (B) UMAP view of 40 065 single cells, colour-coded by assigned cell type. (C) Marker gene expression for each cell type, where dot size and colour represent the percentage of marker gene expression (pct. exp) and the averaged scaled expression (avg. exp. scale) value, respectively. (D) Composition and proportion of cells in AS and control groups. IL, interspinous ligament; LF, ligamentum flavum; SL, supraspinous ligament; SP, spinous process; UMAP, uniform manifold approximation and projection.

We compared the proportion of each cell cluster in the different sample sets and the total population of cells (figure 1D) and observed that the average proportion of stromal cells in the AS group was much lower than that in the control group; this was probably related to the decrease in cell number caused by the strong ossification activity of AS. Neutrophils constituted the majority of all mixed cells.

ScRNA-seq reveals transcriptional features of stromal subpopulations

According to the method described by Baccin et al,14 we divided all stromal cells into three types containing a total of eight subgroups, two Cxcl12-abundant-reticular (CAR) cell types, expressing adipocyte (Adipo-CAR) or osteolineage genes (Osteo-CAR), and MSCs (figure 2A, online supplemental figure S2A). By analysing the marker genes of each subgroup, we found that Osteo-CAR-C5 (Osteo-CAR-Cluster5) exclusively expressed TNFRSF11B and MMP3 (figure 2B), known osteogenic differentiation and cartilage degradation factors.15 16 MSC subpopulations expressed the genes THY117 and NT5E.18 Moreover, MSC-C1 characteristically existed only in the control group and overlapped the expression of the osteogenic-related characteristic genes of MSC-C2 and Osteo-CAR, while retaining the expression of the stem gene NT5E, which we believe is the transition state of MSC-C2 to Osteo-CAR differentiation. We speculated that MSC osteogenic differentiation in AS is more mature and rapid and does not remain long in the intermediate state. Additionally, Adipo-CAR cells characteristically expressed APOE, ADIPOQ and LPL, which are the main features of adipose tissues and cells. Moreover, they highly expressed IGF2, which promotes tumour cell migration19 and chemotactic cytokines such as CXCL12, which promote the migration of immune cells,20 indicating a strong chemotactic function. Of these, Osteo-CAR-C5 (66%), adipo-CAR (69%) and MSC-C2 (65%) were proportionally upregulated in AS (figure 2C).

Supplemental material

Figure 2

Distinct subclusters of the stroma. (A) UMAP of 368 stromal cells (as in figure 1B), annotated and coloured by their clustering. (B) Marker gene expression for each cell type, where dot size and colour represent the percentage of marker gene expression (pct. exp) and the averaged scaled expression (avg. exp. scale) value, respectively. (C) Proportion of each subcluster of stroma between patients with AS and healthy controls. (D) Upper panel: RNA velocity analysis of MSCs, Osteo-CAR cells and Adipo-CAR cells with the velocity field projected onto the UMAP plot of stroma subclusters from (A). Arrows show the local average velocity evaluated on a regular grid and indicate the extrapolated feature states of the cells. Bottom panel: monocle pseudotime analysis revealing the progression of MSCs, Osteo-CAR cells and Adipo-CAR cells. Trajectory reconstruction of all single cells revealing three branches: Pre-branch, Fate-1 and Fate-2. (E) Heatmap showing the scaled expression of differentially expressed transcription factors in three branches as in (D) catalogued into three major gene clusters (labels on left) that vary as a function of pseudo-time, highlighting specific representative transcription factors in each gene cluster along the right margin. From the centre to the left of the heatmap, the kinetic curve from the Pre-branch along the trajectory to the Fate-1 branch. From the centre to the right, the curve from Pre-branch to the Fate-2 branch. (F) GO analysis of differentially expressed genes associated with the three gene clusters in (E) identified unique response pathways for each branch. (G) The expression of c-JUN, c-FOS, CAVIN1, FOSB and EGR1 in AS and control MSC samples by RT-qPCR. (H) H&E staining and the expression of c-JUN (red), c-FOS (green), and CAVIN1 (yellow) in AS and control spinal enthesis samples by IF. Graph showed mean fluorescence intensity (right). The experiments were repeated in triplicate. The representative results are shown. Scale bar: 200 µm. (I) RNA expression of chemokines for each cell type, where dot size and colour represent the percentage of marker gene (pct. exp) and the averaged scaled expression (avg. exp. scale) value, respectively. *p<0.05; **p<0.01. AS, ankylosing spondylitis; CAR, cxcl12-abundant-reticular; EN, enthesis; GO, gene ontology; IF, immunofluorescence; IL, interspinous ligament; MSC, mesenchymal stem cells; RT-qPCR, real-time quantitative PCR; SP, spinous process; UMAP, uniform manifold approximation and projection.

We used scVelo to determine the transcriptional fate of Osteo-CAR and Adipo-CAR lineage cells (figure 2D, up). Monocle pseudotime analysis further corroborated the MSC-C2 to Osteo-CAR and Adipo-CAR transitions. Notably, the pseudotime trajectory began with MSC-C2 and split into two main branches, with Osteo-CAR and Adipo-CAR placed at opposite divergent ends as two terminally differentiated cell types, which was in line with report of Moustapha Kassem’s21 (figure 2D, down, online supplemental figure S2B, C). To elucidate the molecular dynamics that distinguish these two branches, we analysed gene expression dynamics and found that, along the Fate-1 branch, cluster 1 genes activated at the end of the transition were predominantly involved in the GO terms “bone morphogenesis,” “positive regulation of osteoblast differentiation” and “bone mineralization” (EGR1, FOS, JUN, FOSB, and CAVIN1) which are consistent with the features of osteoblastic differentiation. The Fate-2 branch expressed higher levels of cluster 2 genes enriched for the GO terms “lipoprotein transport,” “response to fatty acids,” “lipoprotein metabolic process,” (CITED4, EPAS1, and PPARG) (figure 2E,F, online supplemental figure S2D). Therefore, these unique transcription factor expression patterns characterise a successful Osteo-CAR differentiation path and a functional difference in the Adipo-CAR subcluster.

The experiments in situ were used to verify the difference in proportions of the key cell subgroups above between AS and normal controls. First, we confirmed that normal control tissues were non-inflammatory by immunohistochemical staining of inflammatory cytokines (IL-6, IL-1β) (online supplemental figure S3A). Immunofluorescence staining was used to verify the difference in density (counts/mm2) of the three key subgroups Osteo-CAR-C5 (15.68 vs 8), adipo-CAR (24.48 vs 15.23) and MSC-C2 (16.33 vs 7.13) between AS and normal controls (online supplemental figure S3B–E).

Supplemental material

C-JUN, C-FOS and CAVIN1 are differentially expressed ex vivo and in situ in patients with AS and promote osteogenesis in vitro

To find relevant transcription factors that may lead to hyperosteogenesis of AS, we conducted an association search for genes with the highest significance ranking in cluster 1 in figure 2E and osteogenic differentiation. Among them, c-JUN,22 c-FOS,23 CAVIN1,24 FOSB25 26 and EGR127 28 are reportedly related to osteogenic differentiation in MSCs. Therefore, we examined the expression of c-JUN, c-FOS, CAVIN1, FOSB and EGR1 in AS and control MSC samples using RT-qPCR and found that c-JUN, c-FOS and CAVIN1, but not FOSB or EGR1, were highly expressed in AS (p<0.05) (figure 2G, online supplemental table S5). IF analysis of the spinal enthesis showed high expression of c-JUN, c-FOS and CAVIN1 in AS, consistent with our RT-qPCR results. Using IF staining, we found that the expression of these indicators was mainly concentrated at the ligamentous end of the attachment point (figure 2H). To validate the function of transcript factors, the MSCs from bone marrow were isolated and identified successfully29 (online supplemental figure S4A, B). The shRNA of the three abovementioned transcript factors was designed and transfected into MSCs (online supplemental table S6). And the results indicated that knocking down of c-JUN, c-FOS and CAVIN1 restrained the ossification of MSCs in vitro (online supplemental figure S4C).

Supplemental material

Supplemental material

Supplemental material

Adipo-CAR is an essential subpopulation related to AS progression

In view of the outstanding recruitment capacity of the Adipo-CAR subpopulation, we analysed the expression of chemokines in each stromal cell subpopulation. We observed a very strong neutrophil recruitment capacity of the Adipo-CAR subpopulation, as it specifically expressed CXCL1/2/3/8 (figure 2I). Therefore, we focused on analysing quantitative and functional differences in neutrophils.

Novel CD99_G1 neutrophil subtype expresses osteogenetic genes in patients with AS

ScRNA-seq analysis indicated that neutrophils were the most prominent cell population (figure 1B). The average abundance of neutrophils was compared with bulk profiles on sorted cells and was distinctly clustered into seven subclusters, including CD99_G1 (CD99+MPO+ELANE+), G1 (MPO+ELANE+), G2 (MKI67+TOP2A+), G3 (MMP8+CHIT1+), G4 (MMP8+S100A9+), G5 (CXCL8+CD53+) and G5_CD16_CD62L (SELL+FCGR3B+) (figure 3A,B, online supplemental figure S5A). In addition to the typical G1 markers, the CD99_G1 subcluster specifically and highly expressed PRTN3, IGFBP7, PEBP1 and BEX3. Among these, PRTN3 is associated with neutrophil function and encodes neutrophil azurophilic granules, whose mRNA expression is correlated with normal-density neutrophil activation30 31 (figure 3B). In contrast, CD99_G1 strongly expressed genes related to migration of neutrophils (SELL,32 figure 3B), suggesting that CD99_G1 might tend to migrate to the ligament attachment site.

Supplemental material

Figure 3

Identification and characteristics of neutrophil subclusters. (A) UMAP of 21 420 neutrophils (as in figure 1B), annotated and coloured by their clustering. (B) Marker gene expression for each cell type, where dot size and colour represent the percentage of marker gene expression (pct. exp) and the averaged scaled expression (avg. exp. scale) value, respectively. (C) Proportion of each subcluster of neutrophils between patients with AS and healthy controls. (D) Application of CytoTRACE to predict the trajectory of the differentiation state of neutrophil subclusters. The predicted differentiation state was visualised in UMAP plots based on the differentiation scores. (E) The differential gene expression in CD99_G1 and G1 subclusters. Purple dots: p<0.05, grey dots: p≧0.05. (F) Violin plots showing the distinct expression patterns of the secretory protein genes in each subcluster. (G) PPI networks of hub transcript factors and target genes in CD99_G1 subclusters. (H) GO and KEGG pathway analyses of significant differential genes of CD99_G1 subclusters between AS and control groups. (I) Dot plots of flow cytometry showing percentage of CD99_G1 neutrophils in the AS and control groups. The experiments were repeated in triplicate. The representative results are shown. (J) Statistic results of flow cytometry experiments. (K) Alizarin red staining reflecting MSCs mineralisation under osteogenic induction conditions with different concentrations of IGFBP2/7 intervention. The experiments were repeated in triplicate. The representative results are shown. *p<0.05; **p<0.01. AS, ankylosing spondylitis; Con, control; GO, gene otology; KEGG, kyoto encyclopaedia of genes and genomes; OS, osteogenesis; PPI, protein–protein interaction; UMAP, uniform manifold approximation and projection.

We compared the differences between various subpopulations of neutrophils in the AS and control groups and found two subclusters, CD99_G1 and G2, that were proportionally elevated in AS, suggesting a possible role in its pathogenesis (figure 3C, online supplemental figure S5B). The stemness results for cell subpopulations in figure 3D and online supplemental figure S5C suggest a developmental maturation trend for these cell subpopulations. Since the largest difference in the proportion of the CD99_G1 subcluster was observed between the AS and control groups, we further analysed the differential genes of CD99_G1 and G1 and found that PRSS57, IGFBP7, PRTN3, MANF and SPINK2 were highly expressed in CD99_G1 (figure 3E).

To further explore the potential regulatory role of CD99_G1 neutrophils and other cells in AS, we analysed the mRNA expression of secreted proteins (figure 3F). Among the top-ranked exocrine proteins that were specifically highly expressed in the CD99_G1 subgroup, a significant number are reportedly associated with osteogenic differentiation. To determine the internal mechanisms underlying the elevation of osteogenic cytokines in CD99_G1 neutrophils in AS, SCENIC analysis was applied (online supplemental figure S5D). We found that HOXA10 and SOX4 were strongly expressed in subcluster CD99_G1 (figure 3G, online supplemental figure S5E). In particular, HOXA10 is well related with the osteogenic regulatory gene IGFBP7 (figure 3G). We confirmed the coexistence of CD99_G1 neutrophils with exocrine proteins IGFBP7 and IGFBP2 by IF (online supplemental figure S5F).

We sought to further explore the functions of neutrophils by enriching for the differential genes between subcluster CD99_G1 of the AS and control groups. GO analysis indicated enrichment of neutrophil degranulation, ATP synthesis-coupled proton transport, ATP metabolic process and IRE1-mediated unfolded protein response (figure 3H). KEGG pathway analysis suggested that differential genes were enriched in oxidative phosphorylation, phagosomes and neutrophil extracellular trap formation (figure 3H), which are linked to inflammation. Flow cytometry results confirmed that the percentage of CD99_G1 neutrophils was elevated in the AS group with an average of 52%, compared with 22% in the healthy controls (p<0.01) (figure 3I,J). Furthermore, stimulation tests indicated that these two secreted proteins could enhance the osteogenesis of MSCs (figure 3K).

Monocyte dendritic cell progenitor-recombinant adiponectin receptor 2 (MDP-ADIPOR2) monocytes contribute to inflammation in AS

We referred to Hettinger et al,33 Reyes et al,34 and Jakubzick et al 35 for the definition of marker genes in monocytic lines and classified them into seven subclusters, including GMP (CD34+GATA2+), MDP (MS4A3+ELANE+MPO+), CCR2-CX3CR1-Monocyte (CX3CR1+), DC (CD1C+), MDP-ADIPOR2 (OLR1+MPO+), OLR1-inflammatory-monocyte (OLR1+IL1β+IL6+), and GPR183_DC (GPR183+CD1C+) (figure 4A,B, online supplemental figure S6A, B). We did not find any significant differences between AS and control in the number of monocytic cells in each subpopulation (online supplemental figure S6A, C). Therefore, we attempted to explore the subpopulations of cells associated with AS pathogenesis by analysing the differential gene and functional enrichment between subpopulations (online supplemental table S7). GO and KEGG pathway entries for differential genes in the AS and control MDP-ADIPOR2 subclusters contained multiple pathological processes associated with AS, such as mineral absorption,36 37 IL-17 signalling pathway,38 NF-kappa B signalling pathway,39 tumour necrosis factor signalling pathway40 and cellular response to unfolded protein41 (figure 4C–E). Almost all terms were strongly associated with inflammation. To visually analyse the different GO terms, their related genes, and the relationships between them, we used the ClueGO plug-in of cytoscape (p<0.01). Multiple genes were present in only one AS-related pathological process. For example, PTGS2, FOS, CXCL2 and CXCL8 were associated with genes belonging to two or more terms that might be able to suppress AS inflammation (figure 4F).

Supplemental material

Supplemental material

Figure 4

The heterogeneity of monocytic cells and the proinflammatory phenotype of monocyte dendritic cell progenitor-recombinant adiponectin receptor 2 (MDP-ADIPOR2) subsets. (A) Uniform manifold approximation and projection (UMAP) of monocytic cells (as in figure 1B), annotated and coloured by the clustering. (B) Heatmap of marker gene expression for each cell type. (C) The differential gene expression of MDP-ADIPOR2 cells between the AS and control groups. Purple dots: p<0.05, grey dots: p≧0.05. KEGG pathway (D) and GO (E) analysis of significant differential genes of MDP-ADIPOR2 subclusters between the AS and control groups. (F) ClueGO results of gene enrichment in both the GO and KEGG pathway analyses. AS, ankylosing spondylitis; GO, gene ontology; KEGG, kyoto encyclopaedia of genes and genomes.

Adipo-CAR and CD99_G1 neutrophil cell communication in the AS microenvironment

Elucidating the explicit interaction between immune cells and stromal cells within the osteoimmunological microenvironment would shed light on the pathogenesis of AS.42 To explore the potential crosstalk between Adipo-CAR cells and CD99_G1 neutrophils with other cell types in the AS microenvironment, we used CellPhoneDB, an interactive web application that infers cellular interactions according to the ligand–receptor signalling database.43 The circos plot detected broadcast ligands and demonstrated extensive communication with related receptors (figure 5A), which was consistent with online supplemental figure S7A. Notably, CD99_G1 neutrophils expressed a significantly higher number of receptors corresponding to ligands from Adipo-CAR cells and a significantly higher number of ligands corresponding to receptors expressed by Adipo-CAR cells (figure 5B), indicating a close interaction between Adipo-CAR cells and CD99_G1 neutrophils in the pathogenesis of AS. Our data along with IF experiment suggest that Adipo-CAR cells attract CD99_G1 neutrophils via the CXCL12-CXCR4 axes (online supplemental figure S7B). IGF2/IGF1R and IGF2/IDE were also the major communication relationships between Adipo-CAR and CD99_G1 neutrophils (online supplemental figure S7C). It was further confirmed that the cell number of Adipo-CAR cells was positively correlated with CD99_G1 cells based on IF staining results in situ (p<0.05) (online supplemental figure S7D).

Supplemental material

Figure 5

Cell–cell communication among specific cell types in the spinal ligament attachment site microenvironment. (A) Circos plot showing the potential cell interactions between Adipo-CAR cells (left panel), CD99_G1 neutrophils (right panel), and the other cell types in the spinal ligament attachment site predicted by CellPhoneDB. The node size represents the number of interactions. The width of the edge represents the number of significant ligand–receptor pairs between the two cell types. (B) The dot plot generated by CellPhoneDB showing potential ligand–receptor pairs between CD99_G1 neutrophils (as a receptor in the upper panel and as a ligand in the bottom panel) and all detected cellular types. The dots are coloured according to the mean expression of the ligand–receptor pairs between two clusters, and dot sizes are proportional to the value of −log10 (p value). (C) Predicted regulatory network centred on the microenvironment of the spinal ligament attachment site. CAR, cxcl12-abundant reticular. MSC, mesenchymal stem cell.

CD99_G1 neutrophils showed extensive interactions with MSC-C2 cells, Osteo-CAR-C3 cells, Osteo-CAR-C4 cells and G5_CD16_CD62L neutrophils (figure 5A), indicating that CD99_G1 neutrophils could communicate with these cells. Further analysis suggested that CD99_G1 neutrophils might interact with MSC-C2 cells via the VEGFB-NRP1, OSM-LIFR and OSM-OSMR axes. The most intense ligand receptor relationship between CD99_G1 neutrophils and MSC-C2 cells was through the hepatocyte growth factor-mesenchymal epithelial transition factor (HGF-MET) axis (figure 5B). In contrast, MSC-C2 cells showed the most interactions with Adipo-CAR cells, followed by Osteo-CAR-C4, Osteo-CAR-C3, G5_CD16_CD62L, MSC-C1, Osteo-CAR-C1, Osteo-CAR-C2 and Osteo-CAR-C5 cells, and CD99_G1 neutrophils (figure 5A).

Metabolic characteristics of each cell subclusters

By analysing the metabolic functions of various subpopulations of neutrophils (online supplemental figure S8A), we found that there were differences between fructose and mannose metabolism, galactose metabolism, glycolysis/gluconeogenesis and carbon metabolism in the subcluster CD99_G1. It was also evident that HOXA10 had the strongest regulatory relationship with LDHB in figure 3G, the gene encoding the B subunit of lactate dehydrogenase. There were differences in fructose, mannose and galactose metabolism in subclusters G1, G2 and G5. Finally, there were differences in starch and sucrose metabolism in subcluster G4. We then explored the reason for the differences in cellular function from a metabolic perspective and found that the most significant differences in galactose metabolism, glycolysis/gluconeogenesis and carbon metabolism processes were found in MDP-ADIPOR2, which were downregulated in AS (online supplemental figure S8B).

Supplemental material

Only galactose metabolism and carbon metabolism processes of CD99_G1 neutrophils and MDP-ADIPOR2 monocytes shared differential genes between the AS and control groups (online supplemental figure S8C, D, online supplemental tables S8 and S9). Among these metabolism-related genes, the expression of G6PD was demonstrated to be lower in CD99_G1 neutrophils in patients with AS (15%) compared with healthy controls (46%) by flow cytometry (p<0.01) (online supplemental figure S8E, F).

Supplemental material

Supplemental material

Discussion

AS is a variant of spondyloarthritis characterised by chronic inflammation of the axial skeleton. There is a close relationship between inflammation of the enthesis and disease progression in AS.44 Currently, the specific molecular mechanisms underlying AS inflammation and new bone formation have not been elucidated. Therefore, it is important to understand the interplay between immune cells and bone cells. Although some cells and pathogenic factors at the site of spinal enthesis in AS have been studied,2 45 46 the cellular composition of the microenvironment and the interactions between them remain unclear. In this study, we selected the tissue of ligament attachment points of the spine for single-cell sequencing to elucidate the composition of cell subpopulations and their relationships in situ.

Spinal ankylosis resulting from an increase in pathological new bone formation is a typical feature of AS and the most difficult problem to treat.47 As osteoblasts are the most important effector cells for bone formation and abnormal osteoblast function can lead to disturbances in bone remodelling homeostasis,48 it is critical to study the regulatory mechanisms underlying abnormal osteoblast functions. In our study, we were able to isolate Osteo-CAR cells and found that the Osteo-CAR-C5 subpopulation was closely associated with endochondral osteogenesis and had an increased preponderance in AS, strongly suggesting that spinal ossification in AS is characterised by endochondral osteogenesis.46 49 50 Furthermore, pseudotime analysis indicated that Osteo-CAR-C5 cells were present during the differentiation of MSCs to osteoblasts, suggesting that the Osteo-CAR-C5 subpopulation may be a suitable target to inhibit spine ossification in different disease processes of AS. We also attempted to identify the key transcription factors that regulate osteogenic differentiation of MSCs and found that the expression of c-JUN, c-FOS and CAVIN1 was elevated in AS and confirmed that they all have the ability to regulate the osteogenic differentiation of MSCs through in vitro experiments. Xu et al 11 also found that the pathogenic gene NF-κB abnormally regulates FOS, JUN and JUNB in peripheral mononuclear cells and drives AS progression. Together with our work, these results strongly indicate that these genes may be genetically related to AS.

MRI has allowed the early detection of fatty lesions in the cone of the spine in patients with AS. Moreover, fatty lesions of the bone marrow most likely represent repair tissue as the starting point for syndesmophyte formation.51 Maksymowych et al 52 suggested that fatty degeneration/fat metaplasia is an intermediate step between inflammation and new bone formation in the sacroiliac joints in patients with axial spondyloarthritis. Therefore, inhibition of fatty lesions may play an important role in delaying bone lesions in AS. So we further explored the key transcription factors involved in adipose differentiation of MSCs and found that the Adipo-CAR (IGF2+APOE+LPL+CXCL12+) subcluster had a much higher percentage of cells in the AS group than the control. The cells in this subcluster demonstrated strong chemokine secretion and the ability to recruit neutrophils. Additionally, we identified CITED4, EPAS1 and PPARG as genes that were more highly expressed in the Adipo-CAR subcluster, which were all closely involved in regulating adipose production.53–56 According to our results, cluster 3 was closely associated with the multidirectional differentiation function of MSCs and increased MSC-C2 was the most stemmed transcription factor cluster (figure 2E). JunB mediates the proproliferative and osteogenic effects of Smurf1 on MSCs in cluster 3.57 Notch3 signalling was found to play an important role in all three lineages of MSC differentiation.58–60 Briefly, the MSC-C2 subpopulation may simultaneously mediate AS spinal ossification by upregulating the osteogenic differentiation capacity and enhancing lipogenic differentiation, leading to fat infiltration.

To further explore the effect of Adipo-CAR on the recruitment of neutrophils and identify specific neutrophil subclusters, we subdivided the subclusters of neutrophils and compared their numerical differences between the AS and control groups. The proportion of CD99_G1 neutrophils increased in the AS group. Cell interaction analysis also showed a strong interaction between subclusters CD99_G1 and Adipo-CAR. Therefore, we designed the immunofluorescence experiment to detect the spatial relationship between CD99_G1 neutrophils and Adipo-CAR. We found the colocalisation and positive correlation between these two cell subclusters, which supported the above view. Further analysis showed that CD99_G1 neutrophils mainly expressed various secreted proteins contributing to osteogenesis. CLEC11A, also known as osteolectin, is an osteogenic growth factor highly expressed in bone marrow stromal cells when stimulated by parathyroid hormone.61 Another exocrine protein, NUCB2, not only increases bone density in the femur and lumbar vertebrae of deovalised rats, but also increases alkaline phosphatase activity in mouse MC3T3-E1 preosteoblastic cells and promotes mineralization.62 IGFBP7 is both a negative regulator of RANKL-induced osteoclast formation and bone loss and directly regulates MSC osteogenic differentiation via the Wnt/β-catenin signalling pathway.63 64 Flow cytometry also showed an increased proportion of CD99_G1 neutrophils in AS. The function of IGFBP7/2 in promoting osteogenesis of MSCs has been confirmed in our experiment in vitro. In recent years, more attention has been paid to the role of immune cells in bone formation.42 65 Because of the strong association of these proteins with osteogenic differentiation, CD99_G1 may mediate AS spinal ossification via an exocrine pathway. Using SCENIC analysis, we identified a valuable target molecule, HOXA10, which regulates proliferation and expression of secretory proteins in CD99_G1. HOXA10 expression is restricted to the early stages of myeloid differentiation.66 Lebert-Ghali et al 67 also discovered a vital function of the HOXA10 cluster in the maintenance of hematopoietic stem and progenitor cell activity.

Therefore, CD99_G1 can be recognised as a pro-osteogenesis subpopulation and act as a bridge between lipid infiltration and ectopic ossification. HOXA10 could be a promising therapeutic target for the inhibition of CD99_G1 proliferation and ossification.

GO and KEGG pathway entries for monocytic differential genes in the AS and control groups indicated that the subcluster MDP-ADIPOR2 had a significant proinflammatory feature. ClueGO analysis suggested that PTGS2, FOS, CXCL2 and CXCL8 in monocytic cells were highly correlated with AS-related inflammatory processes involved in MDP-ADIPOR2. Coincidentally, FOS is not only involved in the inflammatory response of AS peripheral blood mononuclear cells but also in the osteogenic differentiation of MSCs.11 Prostaglandin-endoperoxide synthase, also known as cyclooxygenase, is the key enzyme in prostaglandin biosynthesis and plays a critical role in regulating inflammatory responses.68 Therefore, it also showed that the characteristics of MDP-ADIPOR2 subcluster were closely related to AS.

Previous studies have identified an important role for cell metabolism in immune cell activation. In the pathogenesis of autoimmune diseases, the complex microenvironment of the body has an important influence on the metabolism of immune cells during their function and differentiation, thus altering their fate.69 70 Therefore, we analysed whether there are common factors that lead to changes in different subpopulations of immune cells from a metabolic perspective. Our results suggested that galactose metabolism, glycolysis/gluconeogenesis, and carbon metabolism processes were downregulated in CD99_G1 and MDP-ADIPOR2 subclusters in AS. Studies have confirmed that metabolic modes of multiple types of immune cells are transformed from mitochondrial aerobic metabolism in the resting stage to glycolysis to meet the rapidly increasing energy demand of rapid cell proliferation, migration and production of efficacious molecules when stimulated.71–73 Furthermore, galactose reportedly causes ageing of macrophages and inhibits neutrophil function.74 75 Only glycolysis/gluconeogenesis and carbon metabolism processes of CD99_G1 neutrophils and MDP-ADIPOR2 monocytes shared differential genes between the AS and control groups. Therefore, we hypothesised that promotion of glycolysis/gluconeogenesis to accelerate the maturation of disease-related immune cells might provide a novel strategy for disease prevention. Especially, we found G6PD was demonstrated to be lower in CD99_G1 neutrophils in patients with AS by flow cytometry. It is well known that G6PD is a key enzyme in nicotinamide adenine dinucleotide phosphate (NAPDH) synthesis in the pentose phosphate pathway. NADPH promotes the formation of ROS and NETS in neutrophil, which is a marker of mature neutrophil.76 It was indicated that the lower expression of G6PD in AS neutrophils resulted in the obstruction of NADPH production, which makes more of them stay in the early phase. Therefore, upregulation of G6PD expression may be helpful to reduce the proportion of CD99_G1 neutrophils in patients with AS.

In conclusion, our results revealed the cell profiles of spinal enthesis in AS and the correlation between immune cells and osteogenesis, which may contribute to a better understanding of the pathogenesis of AS (figure 5C, online supplemental video 1). However, there were a few limitations to our study. The number of B and T cells captured was small, and neither quantitative nor functional differences were found. Additionally, we quantitatively validated the expression of key cellular subpopulations and genes but did not investigate their functions in vivo and in vitro. The role of such functional cells and molecules should be further investigated in the pathogenesis of AS.

Supplementary video

Data availability statement

Data are available upon reasonable request. Not applicable.

Ethics statements

Patient consent for publication

Ethics approval

This study involves human participants and was approved by Shanghai Changhai Hospital Medical Ethics Committee Approval No.: CHEC2017-163. Participants gave informed consent to participate in the study before taking part.

Acknowledgments

We thank NovelBio Co. for the support of bioinformatics analysis with their NovelBrain Cloud Analysis Platform (www.novelbrain.com). We would like to thank Editage (www.editage.cn) for English language editing.

References

Supplementary materials

Footnotes

  • Handling editor Josef S Smolen

  • XF and CW contributed equally.

  • Contributors Design and perform of the experiments, collection and analysis of the data, and draft of the manuscript: XF and CW. Preparation of single-cell suspension: BJ. Collection of samples: YX and SZ. Drawing and typesetting: JQ. Perform of statistic analysis: ZJ and BZ. Design of the study and revision of the manuscript: WT and WX. WX, as the guarantor, is responsible for the overall content. All authors read and approved the final manuscript.

  • Funding This research was supported by the Shanghai Sailing Program (22YF1459400) and the National Natural Science Foundation of China (82172390).

  • Competing interests None declared.

  • Patient and public involvement Patients and/or the public were not involved in the design, or conduct, or reporting, or dissemination plans of this research.

  • 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.