Single-cell RNA Seq reveals cellular landscape-specific characteristics and potential etiologies for adolescent idiopathic scoliosis.

BACKGROUNDS: Abnormal vertebral growth and development have been found in adolescent idiopathic scoliosis (AIS) patients, and the proliferation and differentiation of bone development-related cells play important roles in its pathogenesis. However, a comprehensive single-cell-level differentiation roadmap in AIS has not been achieved.
METHODS: The present study compared the single-cell level cellular landscapes of spinal cancellous bone tissues between AIS patients and healthy subjects using high throughput single-cell RNA sequencing (scRNA-seq), which covers multiple cellular lineages including osteoblast, chondrocyte, osteoclast and related immunocytes. We constructed the differentiation trajectories of bone development-related cell lineages through pseudotime analysis, and the intercellular-communication networks between bone development-related cells and immunocytes were further developed.
RESULTS: A total of 11 distinct cell clusters were identified according to the genome-wide transcriptome profiles. t-Distributed stochastic neighbor embedding (t-SNE) analysis showed that mesenchymal stem cells (MSC) were classified into three subtypes: MSC-LOXL2, MSC-IGFBP5, and MSC-GJA1. Gene ontology (GO) analysis showed that MSC-GJA1 might possess greater osteoblast differentiation potential than the others. MSC-IGFBP5 was the specific MSC subtype observed only in AIS. There were two distinct gene expression clusters: OB-DPT and OB-OLFML2B, and the counts of osteoblasts derived from AIS was significantly less than that of non-AIS subjects. In AIS patients, MSC-IGFBP5 failed to differentiate into osteoblasts and exhibited negative regulation of cell proliferation and enhanced cell death. CPC-PCNA was found to be the specific chondrocyte progenitor cell (CPC) subtype observed only in AIS patients. The cell counts of OC-BIRC3 in AIS were less than those in controls. Pseudotime analysis suggested two possible distinct osteoclast differentiation patterns in AIS and control subjects. Monocytes in AIS mainly differentiated into OC-CRISP3.
CONCLUSIONS: Our single-cell analysis first revealed differences existed in the cellular states between AIS patients and healthy subjects and found the differentiation disruption of specific MSC and CPC clusters in AIS. Cell communication analysis provided the possible pathogenesis of osteoblast and chondrocyte differentiation dysfunction in AIS.

INTRODUCTION

Adolescent idiopathic scoliosis (AIS) is a three‐dimensional spinal deformity characterized by structural lateral curvature, axial rotation, and deformity in the sagittal profile. It is the most common disorder of spinal deformity in adolescents, with a prevalence from 0.47% to 5.2%. , The progression of AIS with age may lead to adverse long‐term health related quality of life (HRQoL) as well as huge economic burden to family and society. Therefore, identification of its etiology and pathophysiologic mechanism is the key to further understanding of this disease, as well as its treatment and prevention. Several theories have been proposed to explain the etiology of AIS, including genes and heredity, environment, hormones, metabolism, biochemistry, neurological factors, and asymmetric growth. However, the results are controversial and none of these theories could explain all the characteristics of AIS. Asymmetrical vertebrae growth and development along with low bone density has been considered as a potential etiologic factor in the pathogenesis of AIS. , Moreover, vertebrae is the target for other potential etiologic factors like hormones, biomechanics and so on. Therefore, we hypothesized that the proliferation and differentiation of bone development‐related cells might be the core of the etiology of AIS, and all these potential etiologic factors might have synergistic influences on the bone formation, consequently leading to AIS. In the past decades, traditional gene test method, using bulk tissue and cell isolates, has offered great insights in AIS. However, these studies were performed based on an aggregated average of the present cells without considering the heterogeneity among different cell clusters, hence may indicate very little about any specific cell type. Therefore, a detailed specific cell sequencing in AIS patients is critical to understanding the pathogenesis of AIS.

Single‐cell RNA sequencing (scRNA‐seq) was first introduced by Tang et al, and it is considered as a powerful tool for describing highly heterogeneous cell populations. Charles Chan identified human skeletal stem cell (hSSC) using single‐cell RNA sequencing of the growth plate and diaphysis in human fetal femur and confirmed evolutionarily conserved and distinct activated pathways driving SSC‐mediated bone development between mouse skeletal stem cells (mSSCs) and hSSCs. ScRNA‐seq exhibits great advantages of separation of cell type and analysis of the differentiation trajectory in various physiological and pathological development process, , , which provides new insights into the underlying mechanisms of the etiology and progression of AIS. However, a comprehensive single‐cell‐level differentiation landscape in both normal people and AIS patients has not been achieved. Herein, we performed a high throughput scRNA‐seq in control subjects' and AIS patients' cancellous bone samples from vertebrae, and presented a cellular‐state landscape for bone development‐related cells' early differentiation which covered multiple cellular lineages, including mesenchymal stem cells (MSC), osteoblasts (OB), chondrocyte progenitor cells (CPC), chondrocytes (CC), osteoclasts (OC), and others. Then, analysis of pseudotime dynamics was conducted to identify the developmental trajectories of these cell linages and reveal the gene expression dynamics in the process of cell differentiation. Furthermore, we compared the cellular landscape between control subjects and AIS patients to reveal the distinct landscape specific characteristics and potential targets for the underlying etiology and development of AIS. Based on the expression of ligands or complementary receptors we calculated the number of interactions among each cell type and a network was constructed to reveal potential cell‐cell interactions in AIS patients. We believe that elucidating the lineage map and cell‐cell interactions of bone development‐related cells in AIS patients would help enable the molecular diagnosis and treatment of AIS.

METHODS

Tissue sample dissociation

Spinal cancellous bone tissues were obtained from three patients with AIS (all female, aged from 12 to 15) and from other three non‐AIS spinal trauma patients (all female, aged from 10 to 14) who required surgical correction. The detail information of experimental groups and control subjects were shown in Table S1, Supporting Information. Samples were collected with appropriate approval from ethics committee of our university (Local Ethics Committee of Changhai Hospital, SMMU, No. CHEC2017‐163) and followed the detailed procedures below. All the participants gave informed consent in accordance with the 1964 Helsinki declaration.

Briefly, the samples were surgically removed and kept in MACS Tissue Storage Solution (Miltenyi Biotec) until processing. Then tissues were flushed with phosphate‐buffered saline (PBS), which removed any residual blood and soft‐tissue contaminants, and then minced into small pieces (approximately 1‐3 mm) on ice. Single‐cell isolation from the samples was performed with an enzymatic digestion procedure by DNaseII (Worthington) and collagenase IV (Worthington) for 30 minutes at 37°C accompanied with agitation. Following the digestion process, the dissociated single cells were sieved through a 70‐μm cell strainer and centrifuged at 300g for 5 minutes. After removing the supernatant, the pelleted cells were suspended in red blood cell lysis buffer and incubated for 2 minutes to lyse red blood cells. The pelleted cells were washed with RPMI‐1640 medium, followed by resuspension in RPMI‐1640 medium and refiltered through a 35‐μm cell strainer. Viability for dissociated single cells by Calcein‐AM (Thermo Fisher Scientific) and Draq7 (BD Biosciences) staining varied from 65% to 80% for cells profiled, which was satisfactory for the requirements of single cell RNA sequencing. Due to the small number of cells in cancellous bone tissue, previous methods for extracting single‐cell suspension from soft tissue could not obtain enough cells, so we modified the previous protocols and performed it for the first time in cancellous bone tissue. The workflows are shown in Figure S1A,B.

ScRNA‐seq library preparation

The BD Rhapsody system was then used to capture the transcriptomic information of the two single‐cell suspension pools. The pelleted cells were randomly distributed to more than 200 000 microwells through a limited dilution approach, by which a single cell could pair with the bead with oligonucleotide barcodes added into the solution. After the lysis of cell membranes, poly‐adenylated RNA molecules bound to beads and the compounds were collected in a single tube ready for cDNA synthesis. In the reverse transcription process, the 5′ end of each cDNA (the 3′ end of the mRNA transcript) was labeled with a unique molecular identifier (UMI) and cell tag for identifying cell origin. After 18 cycles of preamplification of the adaptor‐ligated cDNA products, amplified cDNA was prepared for next‐generation sequencing. ScRNA‐seq libraries were then harvested from transcripts linked with the 3′ end cell tag amplified in random priming PCR. The library for each sample was sequenced by HiSeq X (Illumina, San Diego, CA) on a 150 bp paired‐end run. BD Resolve analysis pipeline (FASTQ files) were applied to filter the adaptor sequences and remove the low‐quality reads in sequencing libraries. Based on the filtered clean FASTQ data, sequence barcodes with 3′ end UMIs were identified and extracted. We identified 2226 cells from AIS patients and 2241 cells from control subjects, with an average number of reads around 3150 and 2400, respectively, and detected an average of approximately 947 and 767 genes per cell, respectively.

Single‐cell RNA‐Seq data analysis

The Seurat package v2.0.1 was used to analyze the scRNA‐seq expression data from scRNA‐seq libraries. Filtered UMI counts were further normalized according to gene numbers and mitochondrion percentage to create normalized process data. In the following dimensionality reducing procedure, 8500 genes with high expression variability were identified, and 11 principal components were separated using t‐distributed stochastic neighbor embedding (t‐SNE) analysis with the resolution parameter set at 1.5 (as shown in Figure 1A‐C). Cell types were further categorized into distinct clusters by marker gene functional annotation through gene ontology analysis (with databases from NCBI, Uniprot, and AmiGO) and pathway analysis (with the KEGG database), in which the significant P‐value was defined by Fisher exact test, and FDR (false discovery rate) was calculated by the BH test. Based on the gene reads and pathway files analyzed above, Monocle 2 was utilized for pseudotime analysis to place cell clusters along a developmental trajectory according to biological cell differentiation process. Differences in transcriptional profiles (gene expression, specific cell clusters, and pathway analysis) between AIS patients and control subjects and were compared to identify the possible pathogenesis of AIS.

FIGURE 1

Eleven distinct cell clusters were revealed in control subjects and AIS patients. A, ScRNA‐seq library was normalized according to a minimum library size of 200 genes and a maximum of 20% mitochondrial transcript proportion. B, 8500 genes with high expression variability were identified, which are colored by library size, with darker colors indicating larger libraries. C, t‐Distributed stochastic neighbor embedding (t‐SNE) projection where cells sharing similar transcriptome profiles are clustered by colors representing unsupervised clustering results. D, Heatmap analysis using specific gene expression profiles of known cell types. The identity of each cluster was assigned through matching the expression profile with established cell‐specific marker gene expression for B lymphocyte (BLC), T lymphocyte (TLC), mesenchymal stem cell (MSC), chondrocyte progenitor cell (CPC), neutrophil (NP), osteoblast (OB), monocyte (MC), chondrocyte (CC), dendritic cell (DC), osteoclast (OC) and hematopoietic stem cell (HSC). E, The expected gene ontology (GO) terms were used to verify the identity of each cluster. F, Cluster map suggesting the assigned identity for each cluster. G, Cell counts analysis of identified cell types in control subjects versus AIS patients. H, Cell counts proportions of control subjects versus AIS patients shown as a percentage of total specific cell types. nGene, number of genes; percent.mito, percentage of mitochondrial genes

RESULTS

ScRNA‐seq identified multiple cell types in control subjects and AIS patients

Vertebral tissues from AIS patients and control subjects were collected, and the procedures of sample dissociation and single‐cell library preparation and analysis were conducted as shown in Figure S1A‐D. After cell and gene normalization (Figure 1A,B), we quantitatively classified detected cells into 11 discrete cell populations (Figure 1C) based on genome‐wide transcriptome profiles. The results of specific gene expression heatmaps were shown in Figure 1D. According to the different expressions of marker genes, Cluster 2 was identified as mesenchymal stem cell (MSC), which was associated with a high expression of ITGB1, VACAM1, THY‐1 (CD90), NT5E (CD73), and ENG (CD105) , , and negatively associated with hematopoietic stem cell gene markers such as CD34 and PTPRC. Cluster 5 was associated with high expression of osteoblast markers such as RUNX2, COLA1, SPP1, and ANO5, , , which was annotated as osteoblast (OB). Cluster 3 was annotated as chondrocyte progenitor cell (CPC) for the high expression of BIRC5, UBE2C, and DHFR. Cluster 7 was annotated as chondrocyte (CC) for the high expression of COL2A1, SOX9, COL10A1, Aggrecan, and COMP. , , Cluster 9 was annotated as osteoclast (OC) for the high expression of NFKB‐1, NFKB‐2, ATP6V0D1, NFATC1, and OSCAR. , In addition, Clusters 0, 1, 4, 6, 8, and 10 were annotated as immune‐related subsets (Figure 2A‐C). Cluster 0 was associated with high expression of B lymphocyte (BLC) markers CD79A, CD74, MBZ1, HLA‐DRB1, and CXCR4. , Cluster 1 was annotated as T lymphocyte (TLC) for the high expression of CCR5, CD3E, CD25, and IL2RA. , Cluster 4 was annotated as neutrophil (NP) with high expression of CTSG, ELANE, PRTN3, AZU1, and MPO. , , Cluster 6 was enriched expression of TLR2, ITGAM, CD163, and PECAM1, which was annotated as monocyte. , Cluster 8 was annotated as dendritic cell (DC) with high expression of CD1C, FCER1A, SIRPA, and ITGAX. , Cluster 10 was annotated as hematopoietic stem cell (HSC) with high expression of CD34, NT5E, VCAM1, MCAM, and SELE. , ,

FIGURE 2

ScRNA‐seq analysis of immune‐related cell group. A, t‐SNE plots of subclusters of immunocytes. We defined two B lymphocyte subclusters (B‐IGHM and B‐JCHAIN), two T lymphocyte subclusters (T‐TRAC and T‐GNLY), two neutrophil subclusters (NP‐MS4A3 and NP‐CA1), one monocyte subcluster (MC‐VSIR), one dentritic cell subcluster (DC‐NR4A3) and one hematopoietic stem cell (HSCTHBD). B, Cell counts analysis of identified immunocyte subtypes in control subjects versus AIS patients, colored by the source donor. C, Violin plots show the expression distributions of specific marker genes across subclusters of immuocytes. Cell types are represented by different colors

As a further validation, the expected gene ontology (GO) terms were also found in gene‐specific expression cell clusters (Figure 1E). For example, genes specifically expressed in MSC (Cluster 2) were significantly enriched for osteoblast differentiation (P = 4.94E‐05), skeletal system development (P = 0.000778), ossification (P = 0.001159) and collagen fibril organization (P = 1.64E‐11). Specific genes of osteoblasts (Cluster 5) were significantly enriched for osteoblast differentiation (P = 1.25E‐06) and skeletal system development (P = 8.34E‐06). Genes specifically expressed in CPC (Cluster 3) were significantly enriched for positive regulation of collagen biosynthetic process (P = 0.005497), cell cycle (P = 8.95E‐25), and cell proliferation (P = 0.000175). Specific genes of chondrocyte clusters (Cluster 7) were significantly enriched for extracellular matrix organization (P = 8.32E‐25), cartilage development (P = 1.47E‐09), and osteoblast differentiation (P = 2.12E‐07). Specific genes of HSC cluster (Cluster 10) were significantly enriched for hemopoiesis (P = 4.22E‐05), blood coagulation (P = 2.16E‐08), and monocyte differentiation (P = 4.05E‐06). Specific genes of monocyte clusters (Cluster 6) were significantly enriched for immune system processes (P = 1.38E‐187) and osteoclast development (P = 0.002819). Genes specifically expressed in osteoclast clusters (Cluster 9) were significantly enriched for negative regulation of bone mineralization (P = 0.024408). Specific genes of neutrophil clusters (Cluster 4) were significantly enriched for defense responses to bacteria (P = 3.15E‐06) and defense responses to Gram‐positive bacteria (P = 0.000138). Specific genes of TLC cluster (Cluster 1) were significantly enriched for regulation of immune responses (P = 2.17E‐13), T‐cell receptor signaling pathways (P = 5.07E‐12), and T‐cell activation (P = 1.84E‐06). Genes specifically expressed in BLC cluster (Cluster 0) were significantly enriched for positive regulation of immune system processes (P = 1.35E‐14), B‐cell receptor signaling pathways (P = 1.358E‐12), and B‐cell activation (P = 2.59E‐08). Specific genes of DC clusters (Cluster 8) were significantly enriched for regulation of defense responses to antigen processing and presentation of peptide or polysaccharide antigen via MHC class II (P = 1.04E‐12), antigen processing and presentation (P = 1.90E‐10), and toll‐like receptor signaling pathways (P = 1.48E‐09). These analyses strongly indicated that our cell‐type assignments were accurate, as shown in Figure 1F.

From the results of cell counts analysis (Figure 1G,H), which immunocytes in AIS were significantly increased than those in control subjects, we could infer that immune microenvironment might underlie potential etiology of AIS. Therefore, to systematically map the trajectory of bone development in AIS, we extracted cell clusters according to their biological function and formed two main cluster groups for further analysis, including immune‐related cell group and bone development‐related cell group.

In the immune‐related cell group, subclusters were identified according to heterogeneous gene patterns within separated cell clusters (Figure 2A,B). We observed that cluster BLC could be divided into two subclusters, including B‐IGHM associated with enriched expression of IGHM, MS4A1, and CD79A, and B‐JCHAIN, with enriched expression of JCHAIN, IGHG3, and XBP1 (Figure 2C). Cluster TLC were divided into two subclusters, including T‐TRAC enriched for TRAC, CCL5, and CD3E, and T‐GNLY enriched for GNLY, MMP2, and CD2 (Figure 2C). Neutrophil clusters (NP) consisted of NP‐MS4A3 and NP‐CA1 subclusters as well. These subclusters of immunocytes were different types of cells, and constituted the immune microenvironment of bone development, which might impact the proliferation and differentiation of bone development‐related cells.

Comparison of cell subtypes and differentiation trajectories of bone development‐related cells between control subjects and AIS patients

We used Seurat to reperform t‐distributed stochastic neighbor embedding (t‐SNE) analysis for the bone‐development‐related cell group, which yielded three main subcluster gathering regions, including osteoblast, chondrocyte, osteoclast, and respective precursors (Figure 3A,B). Pesudotime analysis was carried out further to map the differentiation trajectories of osteoblast, chondrocyte, and osteoclast between control subjects and AIS patients, respectively.

FIGURE 3

ScRNA‐seq analysis of bone development‐related cell group. A, t‐SNE plots of subclusters of bone development‐related cell group. There are three main subcluster gathering regions, first of which is osteoblast‐derived cell lineage, including three mesenchymal stem cell subclusters and two osteoblast subclusters. The second region consists of chondrocyte‐derived cell lineage, including two chondrocyte progenitor cell subclusters and two chondrocyte subclusters. While the third region consists of osteoclast‐derived cell lineage, which includes one monocyte subcluster and two osteoclast subclusters. B, Cell counts analysis of identified bone development‐related cell subtypes in control subjects versus AIS patients, colored by the source donor

Osteoblast‐derived cell lineage

Bone‐marrow‐derived mesenchymal stem cells (MSC) exhibit a potential differentiation switch between osteogenesis and adipogenesis, and osteogenic differentiation potency determines the growth potential of the skeletal system. According to the different marker genes' expression patterns, we found three distinct osteogenesis‐related mesenchymal stem cells in AIS and control subjects, which involved MSC‐LOXL2, MSC‐IGFBP5, and MSC‐GJA1 (Figure 4D,E). In addition to the common genes expressed in MSC populations, MSC‐LOXL2 was specifically enriched with LOXL2, CCN2, POSTN, CAV1, and MARX5. MSC‐IGFBP5 specifically expressed IGFBP5, HLA‐B, AC005912.1, AC0099560.2, and CXXL12, while MSC‐GJA1 was defined according to its expression of OLFML2B, RPS3AP6, C1QTNF1, CADM1, and GAS1 (Figure 4B). The gene‐related GO enrichment analysis suggested the skeletal system development related and osteoblast differentiation related GO terms were significantly increased in MSC‐GJA1 when compared with MSC‐LOXL2 and MSC‐IGFBP5, which suggested that MSC‐GJA1 possess greater osteoblast differentiation potential than the other two MSC subsets, and be located as downstream of the MSC osteoblast differentiation trajectory. Among the entire detected MSC population, MSC‐LOXL2 and MSC‐GJA1 were identified in control subjects, with cell counts of 423 and 182 in controls, but there were no such cells detected in AIS (Figure 4E). By contrast, MSC‐IGFBP5 was the specific MSC subtype observed only in AIS patients and composes almost 82% (181/220) of the cell counts in the AIS MSC population. With respect to the osteoblast subsets, there were two distinct gene expression clusters: OB‐DPT and OB‐OLFML2B, which were specifically enriched with DPT, CREB3L1, CHAD, CPE, and SPP1, and OLFML2B, RPS3AP6, C1QTNF1, CADM1, and GAS1, respectively (Figure 4A‐D). It was noteworthy that the counts of osteoblasts derived from AIS was significantly less than that of control subjects (2:303 in OB‐DPT and 9:259 in OB‐OLFML2B), which might underlie the lower bone density status in AIS patients than in control subjects (Figure 4E). To address the heterogeneity in the osteoblast differentiation trajectory between AIS and control subjects, Monocle pseudotime analysis was conducted to order related clusters through the osteoblast differentiation process (Figure 4F). The results suggested that there were two possible distinct osteoblast differentiation branches. In control subjects, MSC‐LOXL2 may initiate osteogenic differentiation and then might transfer to a transition state of MSC‐GJA1, and finally differentiate into clusters OB‐DPT and OB‐OLFML2B. Whereas in AIS patients, MSC‐IGFBP5 failed to differentiate into osteoblasts and exhibited negative regulation of cell proliferation and enhanced cell death according to the GO enrichment analysis.

FIGURE 4

Subclusters of osteoblast‐derived cell lineage. A, Heatmaps show the differential gene expression pattern of each subcluster from osteoblast‐derived cell lineage. Top 20 differential genes of each subcluster are shown. B, Violin plots show the expression distributions of specific marker genes across MSC and osteoblast subclusters. Cell types are represented by different colors. C, t‐SNE plots of subclusters of osteoblast‐derived cell lineage (within the circled area). We defined three MSC subclusters (MSC‐LOXL2, MSC‐IGFBP5 and MSC‐GJA1) and two osteoblast subclusters (OB‐DPT and OB‐OLFML2B). D, Cell counts analysis of identified MSC and OB subtypes in control subjects versus AIS patients (within the circled area). E, Cell counts proportions of control subjects versus AIS patients shown as a percentage of total specific cell types in MSC and osteoblast subclusters. F, Monocle pseudotime analysis suggests two possible distinct osteoblast differentiation branches, first of which indicates the normal osteoblast differentiation process in control subjects, whereas the second branch indicates the failure of osteoblast differentiation in AIS patients

To the best of our knowledge, the present study produced the first evidence of osteoblast differentiation dysfunction of MSC in AIS on the single‐cell level and the consequent significant decrease of osteoblasts may explain the abnormal skeletal development pattern in AIS patients. This differentiation trend was similar to osteoblast development in vivo, indicating the differentiation trajectory could stimulate the biological osteogenic process.

Chondrocyte‐derived cell lineage

As mentioned above, Cluster 3 was annotated as CPC, and two subclusters were identified according to gene expression, including CPC‐SLC2A1 and CPC‐PCNA (Figure 5A‐D). CPC‐SLC2A1 might be related to a greater chondrocyte differentiation predisposition because of the enriched expression of SLC2A1, ANK1, HEMGN, ACLS6, and OSBP2. , , Another CPC subcluster was defined as CPC‐PCNA according to its enriched expression of PCNA, CCNE1, SERPINB1, TUBG1, and CENPU (Figure 5B). Those genes have been verified as associated with abnormal mitotic cell cycle transition, DNA damage response, and mismatch repair, which might play an important role in differentiation surveillance from CPCs to chondrocytes. The differences of gene expression in CPC‐SLC2A1 and CPC‐PCNA indicated the differences of cell functions between these subclusters. Based on our findings, we speculated that CPC‐SLC2A1 play a key role in proliferation and differentiation from CPCs to chondrocytes, while CPC‐PCNA might have a role in dysfunction of chondrocyte differentiation because of its enriched expressions of genes. A cell counts analysis was further carried out according to the different sources of samples (Figure 5E). Interestingly, almost all CPC‐SLC2A1 and CPC‐PCNA were identified in AIS patients compared with control subjects, with cell counts of 350 and 42 in AIS and 114 and 0 in controls, respectively. Our results showed that CPC‐PCNA was the specific CPC subtype observed only in AIS patients. Although CPC‐SLC2A1 was found in both control subjects and AIS patients, the cell counts of CPC‐SLC2A1 were significantly greater than those in controls. Two subclusters of chondrocytes were identified in our scRNA‐seq, including CC‐MPP13 and CC‐MT1G with their enriched expression of MMP13, COL10A1, C2orf40, F13A1, and RPL37A, and MT1G, MELTF, STC2, CNMD, and DDX3X, respectively (Figure 5A‐D). Cell counts analysis according to the different sources of samples showed that both CC‐MPP13 and CC‐MT1G existed only in control subjects (cell numbers: 0 vs 254 and 0 vs 189, respectively), and no chondrocyte was identified in AIS patients (Figure 5E). The loss of chondrocytes in AIS patients indicated the dysfunction of endochondral ossification, revealing the potential pathogenesis of AIS.

FIGURE 5

Subclusters of chondrocyte‐derived cell lineage. A, Heatmaps show the differential gene expression pattern of each subcluster from chondrocyte‐derived cell lineage. B, Violin plots show the expression distributions of specific marker genes across CPC and chondrocyte subclusters. Cell types are represented by different colors. C, t‐SNE plots of subclusters of chondrocyte‐derived cell lineage (within the circled area). We defined two CPC subclusters (CPC‐SLC2A1 and CPC‐PCNA) and two chondrocyte subclusters (CC‐MPP13 and CC‐MT1G). D, Cell counts analysis of identified CPC and CC subtypes in control subjects versus AIS patients (within the circled area). E, Cell counts proportions of control subjects versus AIS patients shown as a percentage of total specific cell types in CPC and CC subclusters. F, Monocle pseudotime analysis suggests two trajectories following chondrocyte differentiation, and both CPC‐SLC2A1 and CPC‐PCNA in AIS patients had trouble differentiating into chondrocytes

Pseudotime analysis was further performed to explore the difference of the differentiation trajectory of chondrocytes between control subjects and AIS patients (Figure 5F). We found two trajectories following chondrocyte differentiation: the first branch consisted of two subclusters of chondrocytes, including CC‐MPP13 and CC‐MT1G. However, the second branch consisted only of CC‐MPP13. The different differentiation trajectory of chondrocytes between control subjects and AIS patients was also examined, and the results showed that both CPC‐SLC2A1 and CPC‐PCNA in AIS patients had trouble differentiating into chondrocytes, leading to the consequent loss of chondrocytes in vertebrae. As mentioned above, CPC‐PCNA was enriched for PCNA15, CCNE116, SERPINB117, TUBG1, and CENPU, which were associated with abnormal mitotic cell cycle transition, DNA damage response, and mismatch repair. Therefore, it was easily understood why CPC‐PCNA had differentiation dysfunction into chondrocytes.

Based on our findings, we speculated that the dysfunctions of cellular differentiation from CPCs to chondrocytes in AIS patients lead to the loss of chondrocytes in vertebrae, which might have adverse effects in endochondral ossification. Consequently, the dysfunctions of chondrocyte differentiation and endochondral ossification might play an important role in vertebral growth and pathogenesis of AIS.

Osteoclast‐derived cell lineage

In the osteoclast differentiation lineage, Cluster 6 was annotated as monocytes (MC), the precursors of the osteoclast, with the specific gene expression of SY K, BST2, HLA‐DRA, ERH, and RPSA. According to the different gene expression, two osteoclasts (OC) subclusters were identified, which included OC‐CRISP3 and OC‐BIRC3. OC‐CRISP3 was enriched for expression of CRISP3, CMTM2, BPI, ALOX5AP, and ADGRG3, while OC‐BIRC3 was enriched with BIRC3, CD163, EIF4A3, EREG, and CCL20 (Figure 6A‐D). There were more MCs identified in AIS than in control subjects, with 45 MCs in AIS and eight in controls, respectively (Figure 6D,E). Consequently, AIS claimed the most part of the OC‐CRISP3 population. However, the cell counts of OC‐BIRC3 in AIS were less than those in controls. Pseudotime analysis suggested two possible distinct osteoclast differentiation patterns in AIS and control subjects (Figure 6F). MC in AIS mainly differentiated into OC‐CRISP3, whereas MC in control subjects could differentiate into two osteoclast clusters including OC‐CRISP3 and OC‐BIRC3, which was consistent with the former analysis, suggesting a potential differentiation defect of osteoclast exists in AIS. , Nevertheless, there has been no sufficient evidence to prove the hypothesis until now. In the present study, neither differentiation disruption or significant decrease of specific osteoclast subcluster was found as well, so we did not list osteoclast and corresponding precursor as target potential etiology in AIS.

FIGURE 6

Subclusters of osteoclast‐derived cell lineage. A, Heatmaps show the differential gene expression pattern of each subcluster from osteoclast‐derived cell lineage. B, Violin plots show the expression distributions of specific marker genes across MC and osteoclast subclusters. Cell types are represented by different colors. C, t‐SNE plots of subclusters of osteoclast‐derived cell lineage (within the circled area). We defined one monocyte subcluster (MC‐SYK) and two osteoclast subclusters (OC‐CRISP3 and OC‐BIRC3). D, Cell counts analysis of identified MC and OC subtypes in control subjects versus AIS patients (within the circled area). E, Cell counts proportions of control subjects versus AIS patients shown as a percentage of total specific cell types in MC and OC subclusters. F, Monocle pseudotime analysis suggests MC in AIS mainly differentiates into OC‐CRISP3, whereas MC in control subjects could differentiate into two osteoclast clusters including OC‐CRISP3 and OC‐BIRC3

Finally, we identified all detected cell subtypes shown in Figure 7A,B, and our results showed that there are 11 distinct cell types existing in AIS patients and control subjects. In particularly, the immunocytes in AIS were significantly increased compared with those in control subjects, and the dysfunctions of cellular differentiation of specific MSC and CPC subtypes were detected. Based on these findings, we speculated that cell communications between bone development‐related cells and immunocytes might impact the differentiation of osteoblast and chondrocyte in AIS, and consequently cause adverse effects in endochondral ossification.

FIGURE 7

The whole detected cellular landscape in control subjects and AIS patients. A, t‐SNE projection of all cells, colored and labeled by different cell subtypes. B, t‐SNE projection of all cells, colored by the source donor. Most cell‐type associated clusters are both made up of control subjects and AIS patients, in particular, immunocytes in AIS are significantly increased than those in control subjects

Cell communications between bone development‐related cells and immunocytes

To study the intercellular communication network among detected cell clusters in vertebral tissues, the ligand‐receptor pairs were developed in the form of an interaction heatmap. The results suggested that complex interactions exist among different cell clusters, which could stimulate the differentiation trend of bone development‐related cell types under the influence of complicated cell communications in vivo. The extensive intercellular cross‐talk was shown by the expression pattern of ligand‐receptor pairs in each cell cluster, which included one‐to‐one and one‐to‐many interactions among different cell types, for example, EGFR expressed by MSC‐IGFBP5 could bind only the ligand of TGFβ1 from MC‐SYK, whereas the MSC receptor FR2 could bind ligands of APP, ANXA1, and CAMP from diverse cell populations.

In consideration of the significant increase in the amount of immunocytes, we explored the potential immunological pathogenesis of osteoblast and chondrocyte differentiation disruption in AIS, which were supported by previous studies. , , As a result, cognate receptors for which broadcast ligands were detected within MSC and CPC clusters. The cell communications between MSC and CPC subtypes and immunocytes were shown in Figure 8A,B. while the protein expression heatmap in related cell subtypes was verified in Figure 8C.

FIGURE 8

Cell communications between bone development‐related cells and immunocytes. A, The interaction heatmap suggests the ligand‐receptor pairs existed between MSC, CPC subclusters and immunocytes. B, The illustration of ligand‐receptor interactions within MSC and CPC subclusters. C, The ligands and receptors expression heatmap within MSC, CPC and immunocyte subclusters. D‐F, Comparison analysis of hallmark pathway expression between different MSC and CPC subclusters

When compared with the other two MSC subclusters, MSC‐IGFBP5 showed significant enrichment for LTβR interactions of LTβ‐LTβR with B‐IGHM and TNF‐LTβR with T‐GNLY. Secretion of lymphotoxin β and TNF superfamily cytokines by lymphocytes and the interaction with the lymphotoxin β receptor on stromal cells have previously been described, , which suggested that the gathering of lymphocytes and cytokines and chemokines release mediate the acute inflammatory burst in various diseases. , , The present study now localized this interaction to an MSC‐IGFBP5 subset and lymphocytes, which may lead to immunocytes recruitment and MSC inflammation and underlie the potential immunological pathogenesis of differentiation dysfunction in the AIS MSC population. Meanwhile, we also identified putative interactions between MSC‐IGFBP5 and T‐GNLY that might account for the insufficient immune surveillance of abnormal MSC in AIS. The enrichment expression pattern of HALE in MSC‐IGFBP5 and the HLAE‐NKG2A/ HLAE‐NKG2A:CD94 ligand‐receptor complex suggested a role for T‐GNLY in regulating MSC‐IGFBP5 switch between proliferation and apoptosis. The enriched NKG1A:CD94 on T lymphocytes has been reported to impair IL2 receptor‐induced proliferation and subsequently affect cytotoxicity and cytokine responses when interacting with HLAE. , We found this interaction also existed in the AIS immune microenvironment, probably leading to MSC‐IGFBP5 abnormal proliferation and escape from immune surveillance. In addition, the cognate receptors for growth factor, EGFR and NRP1, were significantly enriched in the MSC‐IGFBP5 subset (Figure 8A) and could interact with VEGFβ expressed by B‐JCHAIN and T‐GNLY and with HBEGF and AREG expressed by dendritic cells, which might explain the excessive proliferation and differentiation dysfunction of MSC‐IGFBP5 in AIS. ,

Analysis of hallmark pathway gene expression demonstrated that the three detected MSCs clusters showed different signal activation pathways. A direct comparison of MSC‐IGFBP5 vs MSC‐LOXL2 and MSC‐GJA1 revealed that the significantly increased pathways in the AIS MSC cluster were the estrogen signaling pathway, NF‐κB signaling pathway, toll‐like receptor signaling pathway, and MAPK signaling pathway, while the TGF‐β signaling pathway was significantly downregulated (Figure 8D,E). Indeed, earlier studies found that the NF‐κB signaling pathway played a key role in suppression of osteoblast differentiation and bone formation in inflammatory bone diseases, , suggesting that NF‐κB and TGF‐β pathways might underlie the osteoblast differentiation dysfunction and low bone density in AIS patients. Taken together, the intercellular cell communications and comparisons in GO and pathway analysis indicated that MSCs of AIS were remodeled to downregulate their osteoblast differentiation through immunocytes recruitment and interaction. Although there was a study reporting immune system abnormalities in AIS paraspinal muscle, , the present study was the first study to report the potential T lymphocyte‐mediated pathogenesis, and further research is necessary to verify this hypothesis.

Identifying the intercellular interactions might explain the chondrocyte differentiation dysfunction in AIS, so we performed the cell communications analysis between CPCs and immunocytes simultaneously. Notably, significant enrichment of IDE‐CCL3, PRPRC‐CD22, NECTIN1‐CD96, and SPN‐SIGLEC were found in CPC‐PCNA. In particular, previous studies revealed that an insulin‐degrading enzyme coded by IDE could terminate insulin activity and participate in intercellular energy management, and decreased expression of IDE could lead to elevated CCL3 levels.

In consideration of the distribution pattern of CPC‐SLC2A1 in AIS and control subjects (350:42), we speculated that CPC‐SLC2A1 could not entirely represent the properties of CPCs in AIS. Therefore, analysis of hallmark pathway gene expression was conducted in CPCs in AIS and those in control subjects, and the results suggested that the HIF signaling pathway, FoxO signaling pathway, and mTOR signaling pathway were significantly increased in CPCs of AIS. Interestingly, the previous study reported enhanced HIF‐1 signaling in chondrocytes resulting in skeletal dysplasia, due to an intracellular energy deficit caused by decreased glucose oxidation. , , Therefore, we inferred that the dysfunction of chondrocyte differentiation in AIS might be attributed to the prolonged HIF‐1α signaling pathway in CPCs, and the consequent oxidative phosphorylation decrease has been verified in our pathway comparison analysis (Figure 8F).

In conclusion, the cell communications between bone development‐related cells and immunocytes implied the potent influence of immunocytes on proliferation and differentiation of specific MSC and CPC subtypes in AIS: (a) The interactions between MSC‐IGFBP and lymphocytes leaded to MSC inflammatory burst, abnormal proliferation and escape from immune surveillance, which might underlie the potential immunological pathogenesis of AIS. (b) Chondrocyte differentiation could be blocked in AIS because of a prolonged HIF‐1α signaling pathway, consequently resulting in the intracellular energy deficit.

DISCUSSION

Altogether, the present study firstly illustrated a comprehensive single cell‐level differentiation roadmap in AIS and revealed cellular‐state landscape differences of bone development‐related cell differentiation between control subjects and AIS patients. We explored a novel cancellous bone tissues dissociation method in preparing single cell samples, which was different from the traditional protocol used in soft tissue bulk, and successfully extracted 2226 cells from AIS patients and 2241 cells from control subjects. Although limited by the properties of cancellous bone, the captured cell was enough to satisfy the requirements of single cell RNA‐seq, verifying the protocol could be effectively employed in bone tissue analysis. Nineteen distinct cell subtypes were identified according to specific gene expression and GO analysis, and the results were consisted with the physiological development in vivo, indicating our study was accurate and might add to our further understanding of the potential targets for AIS etiologies.

Osteoblast‐derived MSC determines the osteogenesis potential in the skeleton system development and plays the pivotal role in various dysplasia diseases. In our study, we found three MSC clusters including MSC‐LOXL2, MSC‐GJA1 and MSC‐IGFBP5, of which the cluster of MSC‐IGFBP5 was unique to AIS. The gene expression of LOXL2, GJA1 and IGFBP5 in MSC have been reported in previous studies. , , In particular, Igfbp5 was used to define MSC populations in periodontium and dental pulp and could enhance osteo/dentinogenic differentiation of MSC in an inflammatory niche. However, our monocle pseudotime analysis revealed that differentiation disruption of MSC‐IGFBP5 in AIS may result in the loss of osteoblasts, and we speculated that might be due to the unique immune environment in AIS. The cell communications between MSC‐IGFBP5 and immunocytes, including LTβ‐LTβR, TNF‐LTβR, and HLAE‐NKG2A, might lead to MSC inflammation burst and insufficient immune surveillance. As a result, the estrogen signaling pathway, NF‐κB signaling pathway, toll‐like receptor signaling pathway, MAPK signaling pathway, and TGF‐β signaling pathway in MSC‐IGFBP5 were abnormal activated in MSC‐IGFBP5, and these may finally underlie the potential immunological pathogenesis of MSC differentiation dysfunction in AIS. These signaling pathway were previously reported to play essential roles in AIS pathogenesis and accompanied osteopenia, , which verified our findings.

Endochondral ossification is the main development model in spine, and hence the differentiation of chondrocyte is of vital importance, but our study found that AIS patients also have defects in this process. Chondrocyte progenitor cell (CPC), which has been identified recently, was divided into two subtypes in our analysis, including CPC‐SLC2A1 and CPC‐PCNA. The gene deletion of SLC2A1 has been proved to severely impair chondrocyte proliferation and hypertrophy, suggesting the potent role of SLC2A1 in cartilage development, while PCNA was proved related with cell proliferation, so both gene‐identified CPC subtypes in our study were consisted with their progenitor roles. Interestingly, we identified more CPCs in AIS than control subjects, whereas chondrocytes were barely detected in patients, which strongly indicated the differentiation dysfunction of chondrocytes in AIS patients, which was accordingly verified by the pseudotime analysis. The cell communication analysis between CPCs and immunocytes suggested significant enrichment of IDE‐CCL3, and insulin‐degrading enzymes coded by IDE could terminate insulin activity and participation in intercellular energy management. The previous study reported that enhanced HIF‐1 signaling in chondrocytes could lead to decreased glucose oxidation and consequently give rise to skeletal dysplasia, which was verified by our analysis of hallmark pathways. Therefore, the dysfunction of chondrocyte differentiation in AIS might be attributed to the prolonged HIF‐1α signaling pathway in CPCs. The intracellular energy deficit induced dysfunction of chondrocyte differentiation and further caused endochondral ossification defects in AIS.

There were some limitations existed in the present study: (a) though we have modified the dissociation protocol, the sample size of captured cells was relative small when compared with soft tissue, which partly was due to less cell counts and special cell distribution pattern in cancellous bone; (b) all the results in our study were based on ingle‐cell RNA sequencing data, lacking of the further experimental validation of proteomics.

CONCLUSION

Our study firstly reported a map of the cellular landscape in AIS using single‐cell RNA sequencing, revealing differences in bone development‐related cell differentiation and cell communications between AIS patients and control subjects. Based on the comparisons of intercellular communications between AIS patients and healthy subjects, several potential pathogenesis of AIS were speculated: (a) The interactions between MSC‐IGFBP and immunocytes in AIS might lead to immunocytes recruitment and MSC inflammation, which could underlie the potential immunological pathogenesis of MSC differentiation dysfunction in AIS. (b) The blocked oxidative phosphorylation in CPC might be due to the abnormally prolonged HIF‐1α signaling pathway, leading to intracellular energy deficit and CPC differentiation dysfunction, and eventually giving rise to defects of endochondral ossification in AIS. Our study might offer new insights into etiology of AIS and potential targets for further treatment.

Figure S1 ScRNA‐seq profiling of cells from the healthy human and AIS vertebra tissues. A, Overview of the single‐cell isolation and analysis workflow. B, Workflows of the dissociation of vertebra tissues and the preparation of single‐cell solutions for scRNA‐seq. C, The generation of scRNA‐seq cDNA expression libraries using BD Rhapsody system. D, Data analysis, including cell clustering and cell counts analysis, pseudotime analysis and cell communication analysis.

Click here for additional data file.

Table S1 Information of AIS patients and control subjects in our study.

Click here for additional data file.