Skip to main content

Transcriptional profiling of human cartilage endplate cells identifies novel genes and cell clusters underlying degenerated and non-degenerated phenotypes

Abstract

Background

Low back pain is a leading cause of disability worldwide and is frequently attributed to intervertebral disc (IVD) degeneration. Though the contributions of the adjacent cartilage endplates (CEP) to IVD degeneration are well documented, the phenotype and functions of the resident CEP cells are critically understudied. To better characterize CEP cell phenotype and possible mechanisms of CEP degeneration, bulk and single-cell RNA sequencing of non-degenerated and degenerated CEP cells were performed.

Methods

Human lumbar CEP cells from degenerated (Thompson grade ≥ 4) and non-degenerated (Thompson grade ≤ 2) discs were expanded for bulk (N=4 non-degenerated, N=4 degenerated) and single-cell (N=1 non-degenerated, N=1 degenerated) RNA sequencing. Genes identified from bulk RNA sequencing were categorized by function and their expression in non-degenerated and degenerated CEP cells were compared. A PubMed literature review was also performed to determine which genes were previously identified and studied in the CEP, IVD, and other cartilaginous tissues. For single-cell RNA sequencing, different cell clusters were resolved using unsupervised clustering and functional annotation. Differential gene expression analysis and Gene Ontology, respectively, were used to compare gene expression and functional enrichment between cell clusters, as well as between non-degenerated and degenerated CEP samples.

Results

Bulk RNA sequencing revealed 38 genes were significantly upregulated and 15 genes were significantly downregulated in degenerated CEP cells relative to non-degenerated cells (|fold change| ≥ 1.5). Of these, only 2 genes were previously studied in CEP cells, and 31 were previously studied in the IVD and other cartilaginous tissues. Single-cell RNA sequencing revealed 11 unique cell clusters, including multiple chondrocyte and progenitor subpopulations with distinct gene expression and functional profiles. Analysis of genes in the bulk RNA sequencing dataset showed that progenitor cell clusters from both samples were enriched in “non-degenerated” genes but not “degenerated” genes. For both bulk- and single-cell analyses, gene expression and pathway enrichment analyses highlighted several pathways that may regulate CEP degeneration, including transcriptional regulation, translational regulation, intracellular transport, and mitochondrial dysfunction.

Conclusions

This thorough analysis using RNA sequencing methods highlighted numerous differences between non-degenerated and degenerated CEP cells, the phenotypic heterogeneity of CEP cells, and several pathways of interest that may be relevant in CEP degeneration.

Introduction

Low back pain is a common and potentially debilitating musculoskeletal disorder. It has been a leading cause of disability worldwide for decades [1] and is predicted to affect up to 80% of individuals in their lifetime [2] and afflict ~20% of older adults annually [3]. Cases relating to back pain are estimated to cost over $100 billion annually in the USA alone due to lost wages and medical expenses [4]. The socioeconomic challenges presented by back pain are projected to increase overtime [5]; thus, it is imperative to elucidate the causes of back pain and to improve therapeutics and preventative measures. Many cases of back pain originate from degenerated intervertebral discs (IVD) [6,7,8]. The IVDs are complex fibrocartilaginous pads that consist of a central gelatinous nucleus pulposus (NP) encapsulated peripherally by concentric collagen rings of the annulus fibrosus (AF) and axially by thin cartilage endplates (CEP). The healthy IVD provides pain-free flexibility to the spine, but it is susceptible to degeneration characterized by a loss of tissue integrity that decreases IVD height, alters IVD mechanics, and facilitates the invasion and sensitization of blood vessels, nerves, and immune cells that generate pain [9,10,11,12,13].

IVD degeneration is associated with a combination of external, lifestyle, and genetic factors [14], with one such factor being injury or degeneration of the CEPs [15]. Structurally, the CEPs encapsulate the NP to separate it from the subchondral bone and maintain a high internal pressure responsible for the disc’s mechanical function [15]. CEP defects can depressurize the NP and are often associated with altered disc mechanics that may initiate IVD degeneration [9, 16, 17]. CEP injury is also a strong predictor of pain [18], possibly by facilitating neurovascular ingrowth [10, 13, 19] and direct interactions between the NP and extradiscal tissue that can cause painful spinal cord sensitization [11] or Modic changes [20, 21]. The CEPs are also the primary route for nutrient transport to sustain the cells of the central region of the disc. Reductions in CEP permeability caused by increased calcification and other matrix compositional changes can disrupt nutrient supply to the central region of the disc and are associated with cell death and increased risk of IVD degeneration [22,23,24]; this could also complicate the efficacy of regenerative therapies for IVD degeneration that target the NP [25]. Therefore, it is evident that CEP degeneration can have a meaningful impact on the progression of IVD degeneration and back pain.

The phenotype and functions of the resident CEP cells and how these change during degeneration are less well characterized. Healthy CEP cells can regulate the tissue extracellular matrix by secretion of collagen II and proteoglycans [26] but upregulate catabolic enzymes during degeneration [27, 28]; loss of CEP cellularity is also associated with matrix degradation in vivo [29]. Degenerated CEP cells also upregulate inflammatory cytokines and may contribute to tissue calcification via the expression of hypertrophic and osteogenic markers [27, 30,31,32,33]. These cells are reminiscent of articular chondrocytes [34], and many studies assess the degenerative status of CEP cells using common markers associated with cartilage health and disease [27, 32, 35]; additionally, several mechanistic studies exploring non-traditional markers in CEP cells are motivated by previous findings in articular cartilage [36, 37]. Despite their similarities, the CEP and articular cartilage are structurally and compositionally distinct tissues [38], and analysis of a subset of cartilage and disc genes revealed phenotypic differences between the two cell types [34]. Taken together, CEP cells should be treated as a distinct cell type that should be better characterized. Additionally, the extent to which CEP cell phenotype shifts during degeneration is not well characterized.

This study aimed to thoroughly characterize the phenotypes of human CEP cells from degenerated and non-degenerated IVD tissue using bulk RNA sequencing followed by single-cell RNA sequencing. First, we used bulk RNA sequencing to compare the broad gene expression of non-degenerated and degenerated human CEP cells and performed a literature search to determine the novelty of the genes from our dataset in CEP cells and other cartilaginous tissues. Next, we used single-cell RNA sequencing to identify and characterize different subpopulations of human CEP cells and to study their expression of the genes in our bulk RNA sequencing dataset. We hypothesized that non-degenerated and degenerated human CEP cells have distinct gene expression profiles and that we will identify numerous genes that were previously unstudied in CEP cells. We also hypothesized that both non-degenerated and degenerated CEP cell samples will consist of numerous cell subpopulations with distinct gene expression and functional profiles. The results of this study may be used to identify novel differences between non-degenerated and degenerated CEP cells and open the door to identifying novel CEP cell-specific markers.

Methods

All reagents were purchased from Thermo Fisher Scientific (Waltham, MA, USA) or VWR International (Radnor, PA, USA) unless otherwise stated.

Cell isolation

Primary cells from de-identified human cadaveric CEP tissue were used for all experiments. Tissue was harvested from cadaveric lumbar spines obtained through Ohio State University’s Cooperative Human Tissue Network under an Institutional Review Board exemption within 48 h of death as described previously [32]. Spines were cleaned of soft tissue and dissected to isolate individual vertebra-disc-vertebra motion segments, which were then sectioned in the sagittal plane and imaged for Thompson grading by two independent investigators. After, the IVDs were separated from the adjacent vertebrae, exposing the CEP tissue. The CEP was carefully scraped away from the vertebrae using a scalpel blade, cut into smaller pieces, and digested sequentially with 0.2% w/v protease from Streptomyces griseus (Cat#P5147; Sigma-Aldrich, St. Louis, MO, USA) for 1 h and 0.2% w/v collagenase II from Clostridium histolyticum (Cat#17101015) for 4 h. The digests were passed through 70-μm strainers and expanded in standard media (high-glucose Dulbecco’s modified Eagle medium (DMEM), 10% fetal bovine serum (FBS), 100 U/mL penicillin/streptomycin, 0.5% fungizone, 50 μg/mL freshly added ascorbic acid) and incubated at 5% CO2, 21% O2, 37 °C for at least one passage. The cells were frozen at − 80 °C in 10% dimethyl sulfoxide in FBS until needed for experiments. CEP cells used in experiments were classified as “non-degenerated” if from an IVD joint of grade ≤ 2 or “degenerated” if from an IVD joint of grade ≥ 4; sample demographics for bulk and single-cell RNA sequencing are provided in Table 1 and Table 2, respectively.

Table 1 Sample demographics for bulk RNA-Seq
Table 2 Sample demographics for single-cell RNA-Seq

Bulk RNA sequencing sample preparation and data analysis

For bulk RNA sequencing, CEP cells (passage ≤ 2) from non-degenerated (N=4) and degenerated (N=4) IVDs were thawed, plated in T75 flasks at a density of 0.5–1 × 106 cells per flask, and expanded until ~80% confluent. During culture, cells were kept at 37 °C, 21% O2, and 5% CO2 and were fed 2–3 times per week with standard media. After expansion, the cells were washed once with 1X PBS and collected in TRIzol reagent (Cat#15596026). RNA was isolated using a PureLink RNA MINI kit (Cat#12183018A), and a preliminary assessment of RNA concentration and purity (260/280 and 260/230 ratios) was performed using a NanoDrop 2000c spectrophotometer. To assess suitability of total RNA for bulk RNA sequencing, RNA Integrity Number (RIN) was measured using BioAnalyzer RNA Nano Kit (Cat#5067-1511; Agilent Technologies, Inc., Santa Clara, CA, USA), and the RNA concentration was fluorometrically quantified using Qubit RNA High Sensitivity Assay (Cat#Q32852).

Samples with RIN >7 were used to generate mRNA sequencing libraries using a NEBNext® Ultra™ II Directional (stranded) RNA Library Prep Kit (Cat#E7760L; New England BioLabs, Inc., Ipswich, MA, USA) and a NEBNext Poly (A) mRNA Magnetic Isolation Module (Cat#E7490; New England BioLabs, Inc.). Briefly, 200 ng of total RNA were used as sample input. RNA fragmentation was set at 10 min, and 12 polymerase chain reaction (PCR) cycles were used in final library generation. Library quantification and characterization were assessed with a BioAnalyzer High Sensitivity DNA Kit (Cat#5067-4626; Agilent Technologies, Inc.) and a Qubit DNA High Sensitivity Assay Kit (Cat#Q32854). Libraries were pooled together with other index-compatible RNA sequencing libraries for sequencing on a NovaSeq 6000 (Illumina, San Diego, CA, USA) paired-end 100bp flow cell to a minimum depth of 20 million clusters per sample.

Basic analysis was performed using a custom in-house pipeline in which raw FASTQ data were aligned to the human reference genome GRCh38 using hisat2 v2.1.0 [39, 40]. Gene-wise counts were generated with featureCounts from the subread package v1.5.1 [41] for genes annotated by ensemble Homo_sapiens.GRCh38.101, counting the primary alignment in the case of multimapped reads. Differential expression analysis was performed using deseq2 [42] and included batch as a covariate to correct for batch biases internally. An rlog normalization of the counts was performed to visualize the results, and a batch correction using removeBatchEffect from limma was performed [43]. A conceptual overview of our bulk RNA sequencing experiment is provided in Fig. 1A.

Fig. 1
figure 1

Conceptual workflow of RNA sequencing experiments. For both bulk RNA sequencing (A) and single-cell RNA sequencing (B), CEP cells from human cadaveric lumbar spines were isolated, expanded, and processed. A list of analyses performed within each dataset are also listed. Figure created with licensed Biorender.com software

Single-cell RNA sequencing sample preparation and data analysis

Passage 1 non-degenerated (N=1) and degenerated (N=1) CEP cells were thawed, plated in T75 flasks at a density of 0.5–1 × 106 cells per flask, and expanded until ~80% confluent. During culture, cells were kept at 37 °C, 21% O2, and 5% CO2 and were fed 2–3 times per week with standard media. After expansion, the cells were processed according to 10X Genomics instructions [44]. Briefly, the cells were put into suspension, filtered using a 70-μm strainer, rinsed three times with 1X PBS containing 0.04% bovine serum albumin (BSA) (Cat#A9418; Sigma-Aldrich), and re-suspended at a concentration of 700–1200 cells/μL in 1X PBS containing 0.04% BSA.

Single-cell suspensions were processed using a Chromium Next GEM Single Cell 3’ Kit v3.1 (PN-1000268, 10X Genomics, Pleastanton, CA) according to manufacturer’s instructions [45]. For each sample, ~16,500 cells were loaded onto a Chromium Next GEM Chip G (PN-2000177, 10X Genomics) to prepare barcoded single-cell gel beads-in-emulsion (GEMs) with a targeted cell recovery of ~10,000 cells. Barcoded cells in GEMs were lysed and cDNA was synthesized and amplified for 11 cycles via PCR. Quality control for cDNA yield and quality was performed using an Agilent 2100 Bioanalyzer at a 1:10 cDNA dilution; samples had final cDNA concentrations between 48 and 54 ng/μL in 40 μL of buffer (1920–2160 ng total cDNA). For each sample, library construction with 25% of the total cDNA content (~475 ng) was started by fragmenting and end-repairing the sample, then amplifying for 10–12 cycles using PCR. Quality control was repeated on the cDNA fragments using a Bioanalyzer High Sensitivity chip (Agilent Technologies, Inc.) prior to sequencing. Samples were sequenced using a NovaSeq 6000 SP flow cell 100-cycle kit (Illumina) in a Read1:i7:i5:Read2 of 28:10:10:90 bp format according to manufacturer’s recommendations [45] to produce 100-bp reads for subsequent analysis. In all, 8035 non-degenerated cells were sequenced at a depth of ~28,000 reads per cell and 9122 degenerated cells were sequenced at a depth of ~24,000 reads per cell, both with a mapping rate of 98%.

Sequenced single cell reads in FASTQ format were aligned to a reference transcriptome using the human genome in 10X Genomics Cell Ranger software (Cell Ranger 6.1.2) [46] based on default parameters with include-introns function. Output files from CellRanger count were analyzed using Seurat 4.0.4 [47]. Cells expressing 200–7000 unique genes and containing < 10% mitochondrial gene content were selected for further analysis. An averaged expression matrix from each cluster was exported from Seurat using the AverageExpression function, and gene expression values were normalized using a z-score calculation function and plotted using pheatmap function [48]. Gene Ontology (GO) functional enrichment analysis was performed using clusterProfiler 4.6.0 toolkit based on default settings [49]. A conceptual overview of our single-cell RNA sequencing experiment is provided in Fig. 1B.

Bulk RNA sequencing literature review

A PubMed literature review was performed to determine which significant genes from the current bulk RNA sequencing dataset were previously identified and studied in different cartilaginous tissues (CEP, IVD, cartilage). The PubMed database was accessed on April 15, 2023, using the National Library of Medicine’s Entrez Programming Utilities following their data extraction guidelines [50]. Searches were conducted by combining relevant search terms for one tissue and one gene at a time. For each gene, the search terms consisted of the gene symbol and the full gene name from genecards.org. For CEP searches, the search terms included “cartilage endplate”, “endplate”, “endplate chondrocyte”, and “cartilaginous”. For IVD searches, the search terms were “IVD”, “intervertebral disc”, “intervertebral disk”, “degenerative disc disease”, “degenerative disk disease”, “nucleus pulposus”, and “annulus fibrosus”. For miscellaneous cartilage searches, the search terms were “chondrocyte”, “cartilage”, “growth plate”, “articular”, “arthritis”, “osteoarthritis”, and “rheumatoid arthritis”. Results were manually reviewed to exclude papers in which the tissue or the gene-of-interest was not the subject of the study; transcriptomic studies, proteomic studies, and review articles that identified or discussed a gene of interest in the tissue of interest were not excluded. The PubMed IDs for all identified papers until April 15, 2023, were saved and the final number of genes with or without prior publications in each tissue of interest were recorded (Table S1).

Results

Bulk RNA sequencing

Numerous genes were differentially expressed between non-degenerated and degenerated CEP samples

We first constructed gene libraries for human CEP cells from non-degenerated and degenerated IVDs using bulk RNA sequencing. 12,965 genes were detected across all samples, of which 76 were significantly differentially expressed (D.E.) between non-degenerated and degenerated CEP cells (q < 0.05). 15 genes were upregulated in non-degenerated CEP cells (fold change ≤ − 1.5), and 38 genes were upregulated in degenerated CEP cells (fold change ≥ + 1.5) (Fig. 2A); 23 additional genes were D.E. but with a fold change magnitude less than 1.5. The fold change and q-values of all significantly D.E. genes are provided in Table S2. The genes in our dataset encoded for proteins that support a broad array of functions, including signaling, molecular transport, transcriptional and translational regulation, metabolism and cell cycling, cytoskeletal and matrix regulation, protein binding and modification, and inflammation. Notably, GALE and HAPLN1, markers involved in extracellular matrix (ECM) synthesis and stabilization [51, 52], were upregulated in non-degenerated CEP cells; conversely, ADAMTS5 and CEMIP, enzymes that break down ECM molecules [28, 53], were upregulated in degenerated CEP cells. This was an expected trend, as excessive matrix breakdown via decreased matrix synthesis and increased MMP and aggrecanase activity is a hallmark of cartilage disease [54, 55].

Fig. 2
figure 2

Bulk RNA sequencing results. A The breakdown of D.E. genes between non-degenerated and degenerated CEP samples is depicted by a volcano plot. Positive F.C. indicates upregulation in degenerated samples, whereas negative F.C. indicates upregulation in non-degenerated samples. Thresholds for upregulation were set at |F.C.| ≥ 1.5 and significance was set at q < 0.05. Boxed numbers indicate number of significantly D.E. genes in each compartment of plot. B The results of the literature review are summarized by a Venn diagram showing the number of genes previously studied in the CEP, the IVD, or cartilage. C Breakdown of the genes previously studied in each tissue of interest. Purple-highlighted genes were previously studied in all three tissues and green-highlighted genes were previously studied in the IVD and cartilage only. F.C., fold change; D.E., differentially expressed

Among the remaining genes, many associated with metabolism and the cell cycle were upregulated in non-degenerated cells, whereas many associated with molecular transport, transcriptional and translational regulation, and cytoskeletal and matrix regulation were upregulated in degenerated cells (Table 3). Specifically, non-degenerated cells were enriched in genes that may be implicated in chondrocyte differentiation and skeletogenesis (FOFX2, PRRX2, S100A2, ID4, KAZALD1), metabolism (ENO1, TST, MPST), matrix production and adhesion (GALE, HAPLN1, LGALS1), and immunogenic functions (CPVL, CYBA). Also enriched in non-degenerated cells were TRPV2, an ion channel that is downregulated in aging cartilage and is necessary for chondrocyte mechanotransduction [56], and AL450405.1, a pseudogene that may be associated with ribosome biogenesis. Several of these genes (MPST, TRPV2) were reported to protect against cartilage calcification in in vivo models of surgically induced osteoarthritis [56, 57].

Table 3 Breakdown of differentially expressed genes with q < 0.05, |fold change| ≥ 1.5 (N = 53 total genes)

Conversely, degenerated cells were enriched in genes associated with signaling and signal transduction (BMPR1B, NTN1, NPR3, ADGRL4, PSD3, SAMD9L, PKD2), endosomal pathways (EEA1, VPS13C, SAMD9L), intracellular transport (CEP290, DYNC2H1, KIAA1109), membrane transport (PKD2, SCN8A, SLC9A7, SLC7A8, AKAP9, ABCC5), cytoskeletal function (TMOD2, DST, AKAP9, FMN1), matrix and adhesion functions (ADAMTS5, CEMIP, THSD4, ITGAV, FAT3, EMB, DSG3), and transcriptional and translational regulation (ATRX, CHD9, CPEB2, GABPB1-AS1, LINC00342, LINC01547, KLF12). Two genes with unknown functions (KIAA1159L, NAALADL2) were also enriched in degenerated cells (Table 3).

Literature search: many differentially expressed genes have not been reported on in CEP cells

CEP cells are significantly understudied compared to their IVD cell counterparts; therefore, we conducted a PubMed literature review to determine whether the genes in our dataset have been previously studied in CEP cells. Of the 76 genes identified, only 2 (ADAMTS5, SULF1) have been studied previously in CEP cells (Fig. 2B). ADAMTS5 is a major aggrecanase in cartilaginous tissues, is upregulated in the degenerated CEP, and is commonly used as a marker of CEP cell degeneration [28, 58]. SULF1 is an enzyme that de-sulfonates heparan sulfated proteoglycans [59, 60], which regulate several signaling pathways such as FGF, Wnt, TGFβ, and hedgehog [60,61,62]; heparan sulfated proteoglycans are also a major constituent of the pericellular matrix proteoglycan perlecan [63]. However, SULF1 was previously studied in a single global in vivo knockout model of disc degeneration [64], and its specific role in CEP cells has not been fully ascertained.

Twenty-nine additional genes were previously studied in cells from the IVD or other cartilaginous tissues, with most being exclusively studied in articular and growth plate chondrocytes (Fig. 2B, Fig. 2C). Interestingly, most genes upregulated in non-degenerated CEP cells were found in this category, and they were mostly associated with metabolic and cell cycling functions (Table S3). Forty-five genes whose relevance in cartilaginous cells has not been previously documented remained (Table S4). Many of these genes were upregulated in degenerated CEP cells and were associated with signaling, molecular transport, transcriptional and translational regulation, and cytoskeletal and matrix regulation. The trends in gene expression between non-degenerated and degenerated CEP cells may point to possible pathways that become dysregulated due to degeneration and may be future targets of investigation.

Single-cell RNA sequencing

Numerous unique cell subpopulations were identified in non-degenerated and degenerated CEP samples

We next performed single-cell RNA sequencing on human CEP cell samples isolated from one non-degenerated and one degenerated IVD joint. Quality control found low enrichment of mitochondrial, ribosomal, and hemoglobin, and platelet RNA, indicating good sample mRNA purity (Fig. S1). Unsupervised clustering identified 11 unique clusters across both samples (Fig. 3A), and differential gene expression analysis indicated that each cluster had a unique gene signature (Fig. 3B, Fig. S2). Clusters were annotated using markers for cell types predicted to reside within the CEP or to be co-isolated during tissue isolation (Table S5, Table S6). We identified 4 distinct chondrocyte clusters, 4 distinct stem cell clusters, an osteoblast cluster, an NP cell cluster, and an AF cell cluster (Fig. 3B). Both samples had a similar compositional breakdown, with chondroprogenitor cells and chondrocytes being the most plentiful cell types while NP cells, AF cells, and osteoblasts were least plentiful (Fig. 3C). Notably, multipotent stem cell and hypertrophic chondrocyte clusters were found exclusively in the degenerated sample. Multipotent stem cells are a class of progenitor cell that can differentiate into several cell types within a particular lineage and includes mesenchymal stem cells [65], though we were unable to resolve the specific cell type within this cluster. Hypertrophic chondrocytes are terminally mature chondrocytes that are upregulated in degenerated cartilaginous tissues and facilitate matrix degradation, matrix calcification, and vascular invasion as part of endochondral ossification [66,67,68].

Fig. 3
figure 3

Single-cell RNA sequencing results. A Unsupervised clustering of non-degenerated (upper) and degenerated (lower) CEP samples reveals sample heterogeneity. Most clusters are shared between both samples, though hypertrophic chondrocyte and multipotent stem cell clusters were only present in the degenerated sample (black circles). Cluster annotations were determined from manually curated markers for different cell types predicted to reside within the CEP. B Heatmap of differentially expressed genes across all clusters. C A breakdown of the cellular composition of each sample, by percentage

Chondrocyte gene and functional enrichment profiles indicate cluster-specific functions

The 4 chondrocyte clusters in our dataset were phenotypically distinct from each other and appeared to have variable functions according to their top differentially expressed genes and GO functional enrichment (Fig. 4, Fig. 5). Whereas hypertrophic chondrocytes are a well-defined cell type, the functions of the other 3 chondrocyte clusters are not as immediately apparent.

Fig. 4
figure 4

Heatmap of differentially expressed genes across all chondrocyte clusters. Cells from the non-degenerated and degenerated samples were pooled in this analysis

Fig. 5
figure 5

Gene ontology comparing functional enrichment of biological processes between different chondrocyte clusters. From left to right: chondrocyte 2 vs. chondrocyte 1, chondrocyte 1 vs. chondrocyte 3, chondrocyte 2 vs. chondrocyte 3. Cells from the non-degenerated and degenerated samples were pooled in this analysis

Cells in the chondrocyte 1 cluster were enriched in pathways such as wound healing, tissue morphogenesis, growth factor stimulus, regulation of cell fate commitment, and cell migration. Wound healing and tissue morphogenesis involve regulation of various cytokine and growth factor pathways to properly coordinate cell migration, proliferation, and differentiation [69]. In support of this, this cluster was enriched in genes that directly or indirectly regulate various signaling pathways (LMO7, CRIM1, FST, IGFBP3, PTX3, EXT1), including the FGF, TGFβ superfamily, hedgehog, and IGF pathways [70,71,72,73,74,75,76]. Based on these observations, this cluster may have a prominent signal regulatory function.

Cells in the chondrocyte 2 cluster expressed genes that regulate stemness (SAT1, EBF1), mitotic activity (CCDC102B, IER2, ID1, ID3), or canonical Wnt signaling (SPON2), which regulates cell proliferation in many cell types [77,78,79]. They also upregulated MEIS2, which regulates chondrogenesis and commonly forms complexes with Pbx1 [80, 81], a patterning factor that is co-expressed with MEIS2 in developing cartilage and is upregulated in proliferating chondrocytes [80, 82]. However, chondrocyte 2 cells were not enriched in markers of G1/S or mitosis (Fig. S3, Fig. S4, Fig. S5, Fig. S6, Fig. S7). Instead, GO analysis showed an enrichment of immune regulation, protein translation, and various metabolic processes. It is possible this is a more quiescent, biosynthetic type of chondrocyte, but the current findings do not support a definitive role for this cluster.

Cells in the chondrocyte 3 cluster expressed genes regulating chondrogenesis and calcification, such as SOX5 [83], ANKH [35, 84], SFRP1 [85], S100A1 [83, 86], and HMGA2 [87]. They also expressed PPFIBP1, which may promote normal long bone development as it is expressed in osteoblasts and has been mapped to a chromosomal locus associated with defects that induce congenital limb shortening [88, 89]. Taken together, cells in this cluster may support a chondrogenic phenotype or may be transitioning to a degenerated pre-hypertrophic phenotype. In support of this, these cells were enriched for pathways related to connective tissue replacement, cytokine production, and positive regulation of chemotaxis, which aligns with the secretion of angiogenic and osteoclastogenic signals by hypertrophic chondrocytes during endochondral ossification to promote vascular and osteoclast invasion into degrading cartilage [68]. Genes associated with the regulation of glycolytic and sugar metabolism pathways were also enriched; this is notable as glycolysis and oxidative phosphorylation dynamics vary at different stages of chondrocyte differentiation [90].

Progenitor cell gene and functional enrichment profiles indicate cluster-specific functions

We next characterized the different progenitor cell clusters in our dataset. Unlike the chondrocytes, visualization of differentially expressed genes showed substantial overlap in the genes expressed by MSCs and proliferating MSCs, and by chondroprogenitors and multipotent stem cells (Fig. 6A). Chondroprogenitor markers did not show enrichment in other progenitor cells but did show enrichment in chondrocyte clusters, supporting that chondroprogenitors fall under a more chondrogenic lineage (Fig. 6B). Overlap in gene expression was observed between MSCs and proliferating MSCs but not in other clusters, supporting that they are closely related cell types (Fig. 6C, Fig. 6D). Multipotent stem cell markers were more uniquely expressed within that cluster but did show some enrichment in degenerated chondroprogenitors and degenerated NP cells (Fig. 6E).

Fig. 6
figure 6

Comparison of progenitor cell clusters in all samples. A Heatmap of differentially expressed genes across all progenitor cell clusters. Cells from the non-degenerated and degenerated samples were pooled in this analysis. BE Dot plots of top differentially expressed markers for B chondroprogenitor cells, C proliferating MSCs, D MSCs, and E multipotent stem cells from the current analysis. Black lines on B separate chondroprogenitors, MSCs, and chondrocyte clusters

Chondroprogenitor cells were enriched in markers for maintenance of tissue ECM (DCN, TIMP1), regulators of oxidative stress (SOD2, SOD3, IGFBP6), heat shock proteins (CLU, CRYAB), and regulators of stemness and proliferation (SPON2, SAT1) (Fig. 6A). CLU, CRYAB, and S100A1 are also upregulated during chondrogenesis [91,92,93] and SAT1 may inhibit inflammation and hypertrophic differentiation of chondrocytes via depletion of polyamines [94]. Furthermore, TIMP1 has an additional antiapoptotic role in various cell types [95]. Taken together, it appeared that chondroprogenitor cells in our dataset regulate tissue homeostasis and cellular stress. GO functional enrichment showed enrichment of lysosomal, endosomal, exosomal, mitochondrial, and immune response pathways (Fig. 7A, Fig. 7B, Fig. 7C). The endosome pathway interacts with lysosome and exosome functions, and the enrichment of lysosomes may suggest an increase in autophagy, which is important for cell survival in the IVD and cartilage [96].

Fig. 7
figure 7

Gene ontology of biological processes comparing different progenitor cell clusters. A MSCs vs. chondroprogenitors. B Proliferating MSCs vs. chondroprogenitors. C Multipotent stem cells vs. chondroprogenitors. D Multipotent stem cells vs. proliferating MSCs. E Multipotent stem cells vs. MSCs. F MSCs vs. proliferating MSCs. For each figure, a description of “Activated” and “Suppressed” is provided

Both MSCs and proliferating MSCs were enriched in pathways relating to cell division, but they may be at different stages of the cell cycle (Fig. 7A, Fig. 7B, Fig. 7D, Fig. 7E). Relative to proliferating MSCs, “normal” MSCs were enriched in pathways such as mitotic cytokinesis, anaphase-promoting complex-dependent catabolic process, chromosome localization, and various pathways related to chromatid/chromosome segregation and separation (Fig. 7F). They upregulated markers associated with microtubule dynamics (STMN1, TUBA1B, PRC1, TPX2), mitotic cyclin regulation (CCNB1, CDKN3), chromosome and chromatid organization (CENPF, TOP2A), and different phases of mitosis (PRC1, STMN1, UBE2S) (Fig. 6D). They also upregulated many genes associated with the different phases of mitosis (Fig. S3, Fig. S4, Fig. S5). Taken together, this may suggest that MSCs are in the process of entering or existing mitosis. Meanwhile, proliferating MSCs were more enriched in pathways such as organic cyclic compound metabolic process, nucleobase-containing compound metabolic process, heterocycle metabolic process, and cellular nitrogen compound metabolic process relative to “normal” MSCs (Fig. 7F); this may suggest an enrichment of nucleotide metabolism, as nucleotides consist of a heterocyclic nitrogenous base, an organic pentose molecule, and a phosphate group. They also upregulated various markers associated with DNA binding, replication, and repair (CLSPN, HELLS, ATAD2, DNMT, RRM1, DEK, NSD2) (Fig. 6C), which are important processes during the DNA synthesis phase (S phase) of the cell cycle [97]; there was a cluster-specific enrichment of G1/S phase markers in proliferating MSCs as well (Fig. S6, Fig. S7). Taken together, this may suggest that proliferating MSCs are more active in the S phase of the cell cycle.

Multipotent stem cells were enriched in numerous ribosomal and mitochondrial pathways (Fig. 7C–E) and upregulated numerous ribosomal markers (RPS5, RPL32, RPL35A, RPS15A, RPL12, RPL29, RPS23, RPL28, RPS12) (Fig. 6E). They may also be involved in further protein processing, as they upregulate SERF2, a marker that promotes protein aggregation, and OST4, which catalyzes the process of co-translational N-glycosylation of proteins to promote proper protein folding and transport [98, 99]. Interestingly, cells in this cluster were also enriched in pathways related to ossification, suggesting a possible role in regulating CEP calcification.

Progenitor cells upregulate genes in bulk RNA sequencing dataset relative to other cell types

Our bulk RNA sequencing results showed that non-degenerated and degenerated CEP cells have different gene expression profiles, though this analysis constitutes the average makeup of the CEP; however, CEP cells are a heterogeneous cell population that all contribute to the overall makeup of those bulk samples. To estimate how the identified subpopulations contribute to overall CEP cell phenotype, we studied their expression of the significantly D.E. genes identified in our bulk RNA sequencing analysis.

Cluster-by-cluster expression of genes upregulated in non-degenerated or degenerated bulk CEP samples were quantified (Table S7) and visualized using dot plots (Fig. 8). We first examined the expression of anabolic (GALE, HAPLN1) and catabolic (ADAMTS5, CEMIP) markers observed in our bulk RNA sequencing results. GALE was enriched in non-degenerated NP cells relative to degenerated NP cells but showed no notable sample-dependent enrichment in the remaining clusters. HAPLN1, however, was upregulated in degenerated chondrocyte 2, chondrocyte 3, proliferating MSC, and chondroprogenitor clusters relative to non-degenerated clusters; these results differed from our bulk RNA sequencing results. In contrast, ADAMTS5 and CEMIP were enriched in hypertrophic chondrocyte and degenerated chondrocyte 1 clusters in cells of the hypertrophic chondrocyte and degenerated chondrocyte 1 clusters, which agreed with our bulk RNA sequencing results; furthermore, ADAMTS5 was enriched in degenerated NP cells relative to non-degenerated NP cells and CEMIP was enriched in degenerated proliferating MSCs relative to non-degenerated cells. There were few notable differences in enrichment between non-degenerated and degenerated clusters for most other genes; CPVL, CYBA, and S100A2 were enriched in most non-degenerated clusters relative to degenerated clusters, whereas BMPR1B, EMB, and PRRX2 were enriched in most degenerated clusters relative to non-degenerated clusters. There was little notable enrichment or downregulation of “non-degenerated” or “degenerative” genes in all chondrocyte clusters. “Degenerative” bulk genes appeared to be downregulated in NP cells but somewhat enriched in osteoblasts relative to “non-degenerated” bulk genes. AF cells, despite accounting for a minor percentage of the cellular composition of our samples, greatly upregulated most “degenerative” bulk genes and downregulated most “non-degenerated” bulk genes.

Fig. 8
figure 8

Dot plot of significantly D.E. markers with |F.C.| ≥ 1.5 identified in bulk RNA-Seq analysis. Vertical line separates the markers on the x-axis between “non-degenerated” bulk genes (left, F.C. < − 1.5) and “degenerative” bulk genes (right, F.C. > + 1.5). Horizontal line separates the clusters on the y-axis between progenitor cells (above) and more differentiated cells (below). D.E., differentially expressed; F.C., fold change

Notably when grouping all progenitor cell clusters together (chondroprogenitors, multipotent stem cells, MSCs, and proliferating MSCs), most “non-degenerated” bulk genes were upregulated and most “degenerative” genes were downregulated. This is notable as mesenchymal stem cells can differentiate and secrete factors with therapeutic effects to augment regeneration of degenerated cartilage and IVD cells [100,101,102,103]. Additionally, a progenitor cell population was previously identified in CEP tissue and may exert pro-chondrogenic and anti-inflammatory effects on other disc cells [104,105,106]. Taken together, this may suggest that a non-degenerated CEP cell phenotype is driven primarily by progenitor cell subpopulations within the tissue.

Degenerated cells were enriched in genes associated with mitochondrial and ribosomal function

After characterizing the heterogeneity of non-degenerated and degenerated CEP cells, we finally studied the effects of degeneration on the different cell subtypes identified in this analysis. GO functional enrichment analysis showed that degenerated chondrocytes and degenerated progenitor cells were enriched in markers for ribosomal, translational, and oxidative phosphorylation-related pathways, whereas non-degenerated clusters were enriched for general cell-cell signaling, communication, and developmental pathways (Fig. 9A, Fig. 9B). Functional enrichment values and gene lists for these pathways are provided in Table S8 (chondrocyte pathways) and Table S9 (progenitor cell pathways). Comparisons between non-degenerated and degenerated cells for individual chondrocyte or progenitor cell clusters showed similar results (Fig. S8-S13). Low ribosomal and mitochondrial RNA content in our quality control analysis suggest these results were not caused by enrichment of ribosomes and mitochondria in our samples (Fig. S1).

Fig. 9
figure 9

Gene ontology of biological processes comparing non-degenerated and degenerated cell clusters. A Non-degenerated chondrocytes 1, 2, and 3 (pooled) vs. degenerated chondrocytes 1, 2, and 3 (pooled). B Non-degenerated chondroprogenitors, MSCs, and proliferating MSCs (pooled) vs. degenerated chondroprogenitors, MSCs, and proliferating MSCs (pooled). For each panel, a description of “Activated” and “Suppressed” is provided

Discussion

Though the CEP has important mechanical and transport functions to regulate IVD degeneration, few studies have characterized the phenotype and functions of the resident CEP cells. While CEP cells resemble articular chondrocytes, several studies have demonstrated tissue and gene-level differences between the two [34, 38]; thus, CEP cells must be more thoroughly characterized to distinguish them from other cartilaginous cell types. Numerous groups have recently used RNA sequencing to characterize the cellularity and degenerative changes within AF and NP tissues [107,108,109,110,111], and Gan et al. [112] has performed single-cell RNA sequencing to examine the cellular profile of CEP tissue from a single donor. To expand on the existing body of knowledge, we performed bulk- and single-cell RNA sequencing on human CEP cells from numerous non-degenerated and degenerated IVDs; the current study investigated broad transcriptomic differences between non-degenerated and degenerated CEP cells, identified novel genes and pathways of interest for future study, and highlighted the phenotypic and functional heterogeneity of CEP cells.

Our bulk RNA sequencing analysis identified numerous differentially expressed genes between non-degenerated and degenerated CEP cells. We identified several anabolic ECM markers (HAPLN1, GALE) that were upregulated in non-degenerated cells and degradative ECM markers (ADAMTS5, CEMIP) that were upregulated in degenerative cells, consistent with prior findings showing an imbalance of ECM marker expression during cartilage and IVD degeneration [27, 54, 55]. Interestingly, all but 2 genes in the current dataset (ADAMTS5, SULF1) were not previously studied in CEP cells, including those that had and had not been previously studied in other cartilaginous tissues (i.e., articular cartilage, growth plate cartilage, IVD). While genes enriched in non-degenerated and degenerated cells both encompassed a broad array of functions, there were disparities in the types of genes enriched in both groups. Many genes enriched in non-degenerated CEP cells were associated with anabolic, metabolic, development, and cell cycling functions (FOXF2, PRRX2, S100A2, ID4, ENO1, TST). FOXF2 and PRRX2 are important factors for craniofacial development, and S100A2 is expressed by chondrocyte-like cells and localizes to sites of calcifying cartilage, suggesting their role in skeletogenesis [111, 113,114,115,116]. ID4 also regulates cell proliferation and differentiation in many developmental processes and is upregulated during chondrogenic differentiation [117, 118]. ENO1 is a major metabolic enzyme that catalyzes the penultimate step of glycolysis and is upregulated by stimulation by the chondrogenic growth factor CCN2 [119]; notably, glycolysis is the primary method of energy production in various cartilaginous tissues due to their hypoxic environments [120, 121]. TST, also known as rhodanese, is another important metabolic enzyme that localizes to the mitochondria and may regulate its function by maintaining redox homeostasis [122]; as mitochondrial dysfunction and reactive oxygen species buildup is observed in osteoarthritic cartilage and aged disc cells [123, 124], this may point to a critical role for TST in regulating metabolism in non-degenerated cells. Several of these genes also have pro-chondrogenic and chondro-protective roles (GALE, HAPLN1, TRPV2, MPST). GALE and HAPLN1 are associated with synthesis and stabilization of proteoglycans within cartilaginous ECM [51, 52]. TRPV2, a mechanosensitive Ca2+ channel, and MPST, a sulfurtransferase that generates the gasotransmitter H2S, are protective against cartilage degradation and calcification in in vivo models of surgically induced osteoarthritis [56, 57]. Interestingly, many of these genes were also previously studied in other cartilaginous tissues.

Degenerated cells were enriched in genes that may regulate chondrocyte hypertrophy and calcification during endochondral ossification, including BMPR1B, FMN1, NPR3, and CEMIP; this is significant as calcification is a major feature of CEP degeneration that decreases tissue permeability and is proposed to occur via endochondral ossification [31, 67]. BMPR1B is a receptor for the chondrogenic growth factor GDF5 and can inhibit hypertrophic differentiation of chondrocytes [125]. FMN1 regulates actin polymerization and regulates chondrocyte hypertrophy and migration during endochondral ossification [126]. NPR3 is a decoy receptor for C-type natriuretic peptide, which is necessary for proper long bone growth [127]; NPR3 deficiency results in elongated bones [127], suggesting a role in regulating endochondral ossification. CEMIP is a degradative enzyme that degrades hyaluronan to facilitate osteoclast and blood vessel invasion into hypertrophic cartilage and is also upregulated in osteoarthritic cartilage [53, 128].

Other genes upregulated in degenerated CEP cells and potentially associated with disease progression and pain were ADAMTS5, NTN1, ITGAV, and THSD4. ADAMTS5 is a major enzyme that degrades the major structural proteoglycan aggrecan and is upregulated in degenerated CEP and other cartilaginous tissues [28]. NTN1 promotes axon guidance in neurons and has been implicated in the ingrowth of sensory neurons in in vivo models of experimental disc degeneration and osteoarthritis [129, 130] and may also promote angiogenesis [131]. Interestingly, nerve ingrowth into the painful disc is usually localized within proteoglycan-depleted regions of tissue [10], highlighting the potential interplay between ADAMTS5 and NTN1. ITGAV is a cell adhesion protein that is expressed by chondrocytes and NP cells [132, 133]. Interestingly, this particular integrin regulates mechanically regulated TGFβ signaling [134] and is upregulated in cartilage and the spine in experimental osteoarthritis and disc degeneration driven by excessive TGFβ signaling [135, 136]. Meanwhile, THSD4 promotes the assembly of microfibrils, which can bind to latent TGFβ binding proteins to regulate the availability of active TGFβ [137, 138]. These results highlight how signaling pathways may become altered in CEP cells in association with degeneration and pain.

Many remaining genes enriched in degenerated CEP cells had not been studied previously in cartilaginous tissues and had signaling, transcriptional and translational regulation, and cytoskeleton and matrix regulation functions. Several of these genes (CEP290, DYNC2H1, PKD2) commonly localize to the primary cilia and are necessary for proper ciliogenesis, intraciliary transport, and ciliary mechanotransduction, indicating a possible significance in signal transduction [139,140,141]. Notably, there appeared to be enrichment of genes associated with endosomal and protein processing pathways (SAMD9L, PSD3, EEA1, VPS13C, KIAA1109, SLC9A7). Endosomes are important for general cell metabolism by functioning as sorting and transport centers for endocytosed and intracellular cargo and preparing it for subsequent processing, degradation, or exosomal re-export [142]. The endosomal pathway is also linked to autophagy, a process of recycling intracellular cargo to promote cell survival and is the common focus in various mechanistic studies of CEP degeneration [105, 143,144,145,146]. Our results appear to corroborate the possible significance of this pathway and highlight that alterations in intracellular signaling and metabolism might be prevalent in degenerated CEP cells.

Interestingly, our bulk RNA sequencing analysis found an enrichment of several cytoskeletal and membrane transport markers with neuron-specific enrichment (SCN8A, FAM155A, KIAA1109, ABCC5, TMOD2) in our degenerated bulk samples. However, we did not identify any neuronal cells despite screening for them in our single-cell analysis. Neural ingrowth into the CEP has been reported in painful IVDs and was found to occur primarily through degenerated, proteoglycan-depleted tissue [10, 13]. Enrichment of these markers in our samples suggests that those markers are naturally expressed by CEP cells, highlighting the use of RNA sequencing to discover novel genes in CEP cells and may suggest the importance of ion channels and cytoskeletal proteins in CEP cell function.

The current single-cell analysis identified various chondrocyte clusters. These were the most plentiful cell types in our samples and each appeared to have distinct functions, including signal regulation, biosynthesis, and pro-chondrogenic or pre-hypertrophic functions. This highlights the cellular heterogeneity within CEP cell samples and presents multiple cellular targets that could contribute to CEP health and degeneration, though it remains to be determined the significance of these different clusters. Notably, Gan et al. [112] also found that the predominant cell types in one CEP sample were 3 primary chondrocyte-like clusters. Among these were “regulatory chondrocytes” that upregulated growth factors and chondrogenic pathway regulators, “homeostatic chondrocytes” that upregulated ECM proteins and circadian rhythm markers, and “effector chondrocytes” that had a more hypertrophic character; this was similar to our findings and supports the presence of these cell subtypes in CEP tissue.

We also identified a chondroprogenitor cluster, two subtypes of mesenchymal stem cells, and another multipotent stem cell cluster. Resident progenitor cells similar to mesenchymal stem cells have been previously identified in articular cartilage [147], the NP [148], and the CEP [104] and are believed to contribute to tissue homeostasis and regeneration; these progenitors have increased chondrogenic potential compared to those isolated from other common stem cell sources and actively proliferate. Mesenchymal stem cells also secrete therapeutic factors that can mitigate tissue injury and inflammation in osteoarthritis and IVD degeneration [100,101,102,103, 105]. In our dataset, our mesenchymal stem cell clusters were actively proliferative whereas our chondroprogenitors upregulated markers associated with ECM maintenance and regulation of oxidative stress and protein aggregation; these results may suggest that progenitor cells responsible for proliferating and secreting therapeutic factors belong to different subpopulations and therefore have different roles in maintaining tissue homeostasis.

Interestingly, our chondroprogenitors were enriched in pathways such as extracellular exosome and those associated with endosomes. Exosomes are small membrane-enclosed structures that can modulate cell activity via the transfer of various biological cargo and are secreted from multivesicular bodies derived from late endosomes [149, 150]. Stem cell-derived EVs have been studied for their therapeutic effects in various diseases, including arthritis and disc degeneration [100,101,102,103]. Notably, exosomes have been isolated from rat CEP cells and were demonstrated to have anti-inflammatory and pro-chondrogenic effects on other disc cells [105, 106], suggesting a possible role for EV-mediated functions for CEP-derived progenitor cells. The enrichment of endosomal-related pathways may also point toward an increase in lysosomal and autophagy-related pathways, thereby contributing to cell survival [96].

Multipotent stem cells were only present in our degenerated sample and were enriched in ribosomal, protein translation, and mitochondrial function pathways. Interestingly, all clusters from our degenerated CEP sample were also enriched in similar types of pathways. Ribosomes are organelles that translate proteins from mRNA transcripts and can therefore regulate protein expression, and their role in regulating cartilaginous disease is gaining new attention [151]. Osteoarthritic and aged chondrocytes exhibit altered ribosome biogenesis and processing, which can regulate their response to oxidative stress and other stimuli [151,152,153]. Whereas chondrocytes typically translate proteins via a 5′ cap-dependent mechanism, chondrocytes exposed to inflammatory stimuli (i.e., cytokines, osteoarthritic synovial fluid) upregulate translation through an alternative IRES-dependent mechanism that is commonly activated in various stress conditions (i.e., endoplasmic reticulum stress, nutrient deprivation, inflammation, apoptosis, mitosis) to produce stress-related proteins [151, 154,155,156,157]. Currently, the association between proper protein translation dynamics and IVD and CEP degeneration is understudied. However, Yang et al. [158] found that several ribosomal proteins and protein translation pathways were differentially expressed between healthy and degenerated IVD cells. This, in coordination with our findings of enrichment of genes associated with protein translation pathways in degenerated CEP cell clusters, indicates ribosome function and protein translation dynamics may be novel factors regulating CEP degeneration.

Mitochondria generate ATP through oxidative phosphorylation and are therefore important for energy metabolism. Articular chondrocytes and NP cells exist in hypoxic environments and preferentially generate ATP via glycolysis [120, 121]. The NP directly underlying the CEP in canine discs has an oxygen tension ~15% that of blood oxygen tension [159], which may suggest the CEP is hypoxic as well and may prefer glycolytic pathways. Though we expanded both non-degenerated and degenerated CEP cells at 21% O2, only degenerated cells were enriched in pathways associated with mitochondrial function and oxidative phosphorylation. Altered mitochondrial function is associated with various diseases, including degeneration of cartilaginous tissues [160, 161]. Some studies have reported that mitochondrial function is reduced in degenerated IVD cells and osteoarthritic cartilage, while upregulation of oxidative phosphorylation has a chondroprotective effect in articular chondrocytes [123, 124]. In contrast, Cisewski et al. reported elevated oxygen consumption rates indicative of increased oxidative phosphorylation in degenerated IVD and CEP cells relative to healthy cells [162]. These reports suggest there could be multiple modes of mitochondrial dysfunction that are associated with IVD degeneration or it may serve as a compensatory mechanism to alleviate degenerative changes.

These results provide a new and substantial body of information expanding our current knowledge of CEP cells, but their impact may be affected by several limitations. We utilized cultured CEP cells instead of freshly isolated cells in these experiments to increase our cell count for RNA sequencing experiments due to low yield of viable cells from freshly isolated human CEP tissue. The results here could be used as a reference for the many experiments that use cultured CEP cells, as monolayer cell culture is a common model used to study CEP cells. However, cell expansion in monolayer could result in a loss of the native CEP cell phenotype through dedifferentiation [163] and may lead to a loss of less adherent cell subpopulations. We attempted to minimize the effects of cell culture by using cells that have been minimally passaged. Though CEP cell phenotype post-expansion was not evaluated prior to RNA sequencing, CEP cells passaged once after isolation continue to express chondrogenic markers ACAN and COL2A1 at similar or elevated levels to articular chondrocytes and NP cells [34], suggesting some maintenance of a native phenotype; furthermore, non-degenerated and degenerated disc cells first expanded in monolayer still have inherent differences post-expansion, suggesting similar differences may be detectable in our dataset [12, 164]. The use of human cells increases the clinical utility of the current work, but it may also increase the variability of the results due to the different demographics of human subjects the samples were sourced from. While this may make it more difficult to identify all differences between non-degenerated and degenerated samples, the future expansion of the current dataset to account for all demographic differences (i.e., sex, age, race, lifestyle) can highlight new findings. We used the macroscopic Thompson scale to grade IVDs and classify CEP cells as “non-degenerated” or “degenerated,” though CEP structure has only a minor influence on scoring and may therefore not reflect the actual degenerative status of our samples [165]. Furthermore, each sample included cells from both the cranial and caudal CEP, which may have different degeneration patterns. Finally, while we did not validate our sequencing results using other quantitative measures such as qRT-PCR or protein assays, the goal of the current work was to explore CEP cell phenotype, highlight novel markers, and identify possible markers that are differentially expressed in non-degenerated and degenerated CEP cells; future experiments will be performed that validate the findings of this sequencing data.

Conclusions

We identified numerous differentially expressed genes between non-degenerated and degenerated CEP cells using bulk RNA sequencing and highlighted many novel genes not previously reported on in the CEP or other cartilaginous tissues. We reviewed the genes in our dataset and pinpointed several pathways of interest, including transcriptional regulation, translational regulation, intracellular transport, and mitochondrial pathways, that may be worth investigating in future experiments. We also characterized the cellular profile of non-degenerated and degenerated CEP cell samples using single-cell RNA sequencing and found numerous subpopulations of chondrocyte-like cells and progenitor cells. Notably, we found that progenitor cells had higher expression of many genes upregulated in non-degenerated CEP cells and lower expression of many genes upregulated by degenerated CEP cells. This work addresses a critical gap in clarifying the phenotype of CEP cells and could act as a stepping-stone to identify novel markers that regulate CEP degeneration, as well as identify novel CEP-specific markers for the development of preclinical models or CEP-targeted therapeutics. This will significantly impact how we understand the role of the CEP in regulating IVD joint degeneration and discogenic back pain.

Availability of data and materials

All figures and tables referenced in this study are available within this published article or available as supplementary files. Raw bulk (GSE242040) and single-cell RNA sequencing data (GSE242443) are publicly available in the GEO database.

Abbreviations

AF:

Annulus fibrosus

BSA:

Bovine serum albumin

CEP:

Cartilage endplate

D.E.:

Differentially expressed

DMEM:

Dulbecco’s modified Eagle medium

ECM:

Extracellular matrix

FBS:

Fetal bovine serum

GEM:

Gel bead-in-emulsion

GO:

Gene Ontology

IVD:

Intervertebral disc

NP:

Nucleus pulposus

PCR:

Polymerase chain reaction

RIN:

RNA integrity number

References

  1. James SL, et al. Global, regional, and national incidence, prevalence, and years lived with disability for 354 diseases and injuries for 195 countries and territories, 1990–2017: a systematic analysis for the Global Burden of Disease Study 2017. Lancet. 2018;392:1789–858.

    Article  Google Scholar 

  2. Andersson G. Epidemiological features of chronic low-back pain; 1999.

    Book  Google Scholar 

  3. Wu A, et al. Global low back pain prevalence and years lived with disability from 1990 to 2017: estimates from the Global Burden of Disease Study 2017. Ann Transl Med. 2020;8:299–9.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Katz JN. Lumbar disc disorders and low-back pain: socioeconomic factors and consequences. J Bone Joint Surg Am. 2006;88:21–4.

    PubMed  Google Scholar 

  5. Ferreira ML, et al. Global, regional, and national burden of low back pain, 1990–2020, its attributable risk factors, and projections to 2050: a systematic analysis of the Global Burden of Disease Study 2021. Lancet Rheumatol. 2023;5:e316–29.

    Article  Google Scholar 

  6. de Schepper E, et al. The association between lumbar disc degeneration and low back pain: the influence of age, gender, and individual radiographic features. Spine (Phila Pa 1976). 2010;35:531–6.

    Article  PubMed  Google Scholar 

  7. Depalma MJ, Ketchum JM, Saullo T. What is the source of chronic low back pain and does age play a role? Pain Med. 2011;12:224–33.

    Article  PubMed  Google Scholar 

  8. van Tulder M, Assendelft W, Koes B, Bouter L. Spinal radiographic findings and nonspecific low back pain: a systematic review of observational studies. Spine (Phila Pa 1976). 1997;22:427–34.

    Article  PubMed  Google Scholar 

  9. Vergroesen PPA, et al. Mechanics and biology in intervertebral disc degeneration: a vicious circle. Osteoarthr Cartil. 2015;23:1057–70.

    Article  Google Scholar 

  10. Lama P, Le Maitre CL, Harding IJ, Dolan P, Adams MA. Nerves and blood vessels in degenerated intervertebral discs are confined to physically disrupted tissue. J Anat. 2018;233:86–97.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Lai A, et al. Spinal cord sensitization and spinal inflammation from an in vivo rat endplate injury associated with painful intervertebral disc degeneration. Int J Mol Sci. 2023;24:3425.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Wiet MG, et al. Mast cell-intervertebral disc cell interactions regulate inflammation, catabolism and angiogenesis in discogenic back pain. Sci Rep. 2017;7:12492.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Freemont AJ, et al. Nerve growth factor expression and innervation of the painful intervertebral disc. J Pathol. 2002;197:286–92.

    Article  PubMed  Google Scholar 

  14. Adams MA, Roughley PJ. What is intervertebral disc degeneration, and what causes it? Spine (Phila Pa 1976). 2006;31:2151–61.

    Article  PubMed  Google Scholar 

  15. Fields AJ, Ballatori A, Liebenberg EC, Lotz JC. Contribution of the endplates to disc degeneration. Curr Mol Biol Rep. 2018;4:151–60.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Wang D, et al. Ex vivo biomechanical evaluation of acute lumbar endplate injury and comparison to annulus fibrosus injury in a rat model. J Mech Behav Biomed Mater. 2022;131:105234.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Adams MA, Freeman BJC, Morrison HP, Nelson IW, Dolan P. Mechanical initiation of intervertebral disc degeneration. Spine (Phila Pa 1976). 2000;25:1625–36.

    Article  PubMed  Google Scholar 

  18. Bailey JF, et al. The relationship between endplate pathology and patient-reported symptoms for chronic low back pain depends on lumbar paraspinal muscle quality. Spine (Phila Pa 1976). 2019;44:1010–7.

    Article  PubMed  Google Scholar 

  19. Stefanakis M, et al. Annulus fissures are mechanically and chemically conducive to the ingrowth of nerves and blood vessels. Spine (Phila Pa 1976). 2012;37:1883–91.

    Article  PubMed  Google Scholar 

  20. Dudli S, Fields AJ, Samartzis D, Karppinen J, Lotz JC. Pathobiology of Modic changes. Eur Spine J. 2016;25:3723–34.

    Article  PubMed  Google Scholar 

  21. Dudli S, et al. Intervertebral disc/bone marrow cross-talk with Modic changes. Eur Spine J. 2017;26:1362–73.

    Article  PubMed  Google Scholar 

  22. Wong J, et al. Nutrient supply and nucleus pulposus cell function: effects of the transport properties of the cartilage endplate and potential implications for intradiscal biologic therapy. Osteoarthr Cartil. 2019;27:956–64.

    Article  Google Scholar 

  23. Roberts S, Urban J, Evans H, Eisenstein S. Transport properties of the human cartilage endplate in relation to its composition and calcification. Spine (Phila Pa 1976). 1996;21:415–20.

    Article  PubMed  Google Scholar 

  24. Bonnheim NB, et al. The contributions of cartilage endplate composition and vertebral bone marrow fat to intervertebral disc degeneration in patients with chronic low back pain. Eur Spine J. 2022;31:1866–72.

    Article  PubMed  Google Scholar 

  25. Huang YC, Urban JPG, Luk KDK. Intervertebral disc regeneration: do nutrients lead the way? Nat Rev Rheumatol. 2014;10:561–6.

    Article  PubMed  Google Scholar 

  26. Antoniou J, et al. The human lumbar endplate: evidence of changes in biosynthesis and denaturation of the extracellular matrix with growth, maturation, aging, and degeneration. Spine (Phila Pa 1976). 1996;21:1153–61.

    Article  PubMed  Google Scholar 

  27. Zhang J, et al. Expression of matrix metalloproteinases, tissue inhibitors of metalloproteinases, and interleukins in vertebral cartilage endplate. Orthop Surg. 2018;10:306–11.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Chen S, et al. Upregulation of tumor necrosis factor α and ADAMTS-5, but not ADAMTS-4, in human intervertebral cartilage endplate with modic changes. Spine (Phila Pa 1976). 2014;39:E817–25.

    Article  PubMed  Google Scholar 

  29. Ariga K, et al. The relationship between apoptosis of endplate chondrocytes and aging and degeneration of the intervertebral disc. Spine (Phila Pa 1976). 2001;26:2414–20.

    Article  PubMed  Google Scholar 

  30. De Luca P, et al. Intervertebral disc and endplate cells response to IL-1β inflammatory cell priming and identification of molecular targets of tissue degeneration. Eur Cell Mater. 2020;39:227–48.

    Article  PubMed  Google Scholar 

  31. Hristova GI, et al. Calcification in human intervertebral disc degeneration and scoliosis. J Orthop Res. 2011;29:1888–95.

    Article  PubMed  Google Scholar 

  32. Lakstins K, et al. Investigating the role of culture conditions on hypertrophic differentiation in human cartilage endplate cells. J Orthop Res. 2020;39:1204–16. https://doi.org/10.1002/jor.24692.

    Article  PubMed  Google Scholar 

  33. Xu H, et al. Intermittent cyclic mechanical tension-induced calcification and downregulation of ankh gene expression of end plate chondrocytes. Spine (Phila Pa 1976). 2012;37:1192–7.

    Article  PubMed  Google Scholar 

  34. Lakstins K, et al. Characterization of the human intervertebral disc cartilage endplate at the molecular, cell, and tissue levels. J Orthop Res. 2021;39:1898–907.

    Article  PubMed  Google Scholar 

  35. Liu MH, et al. Matrix stiffness promotes cartilage endplate chondrocyte calcification in disc degeneration via miR-20a targeting ANKH expression. Sci Rep. 2016;6:25401.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Jiang C, et al. Inhibition of EZH2 ameliorates cartilage endplate degeneration and attenuates the progression of intervertebral disc degeneration via demethylation of Sox-9. EBioMedicine. 2019;48:619–29.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Wang B, et al. miR-142-3p and HMGB1 are negatively regulated in proliferation, apoptosis, migration, and autophagy of cartilage endplate cells. Cartilage. 2021;13:592S–603S.

    Article  PubMed  PubMed Central  Google Scholar 

  38. DeLucca JF, et al. Human cartilage endplate permeability varies with degeneration and intervertebral disc site. J Biomech. 2016;49:550–7.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Gadepalli VS, Ozer HG, Yilmaz AS, Pietrzak M, Webb A. BISR-RNAseq: an efficient and scalable RNAseq analysis workflow with interactive report generation. BMC Bioinformatics. 2019;20:1–7.

    Article  Google Scholar 

  40. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Liao Y, Smyth GK, Shi W. FeatureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.

    Article  PubMed  Google Scholar 

  42. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1–21.

    Article  Google Scholar 

  43. Ritchie ME, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47.

    Article  PubMed  PubMed Central  Google Scholar 

  44. 10X Genomics. Cell preparation for single cell protocols. [Updated 23 Jun 2023]. Accessed from https://www.10xgenomics.com/support/single-cell-gene-expression/documentation/steps/sample-prep/single-cell-protocols-cell-preparation-guide

  45. 10X Genomics. Chromium Next GEM Single Cell 3′ Reagent Kits v3.1 (Dual Index). 2021. www.10xgenomics.com/trademarks.

  46. Zheng GXY, et al. Massively parallel digital transcriptional profiling of single cells. Nat Commun. 2017;8:14049.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Hao Y, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184:3573–3587.e29.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32:2847–9.

    Article  PubMed  Google Scholar 

  49. Wu T, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3).

  50. Sayers E. A General Introduction to the E-utilities. 2009 May 26 [Updated 2022 Nov 17]. In: Entrez Programming Utilities Help [Internet]. Bethesda (MD): National Center for Biotechnology Information (US); 2010. Available from: https://www.ncbi.nlm.nih.gov/books/NBK25497/.

  51. Wen Y, et al. The critical role of UDP-galactose-4-epimerase in osteoarthritis: modulating proteoglycans synthesis of the articular chondrocytes. Biochem Biophys Res Commun. 2014;452:906–11.

    Article  PubMed  Google Scholar 

  52. Urano T, et al. Single-nucleotide polymorphism in the hyaluronan and proteoglycan link protein 1 (HAPLN1) gene is associated with spinal osteophyte formation and disc degeneration in Japanese women. Eur Spine J. 2011;20:572–7.

    Article  PubMed  Google Scholar 

  53. Deroyer C, et al. CEMIP (KIAA1199) induces a fibrosis-like process in osteoarthritic chondrocytes. Cell Death Dis. 2019;10:103.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Le Maitre CL, Hoyland JA, Freemont AJ. Catabolic cytokine expression in degenerate and herniated human intervertebral discs: IL-1β and TNFα expression profile. Arthritis Res Ther. 2007;9:1.

    Article  Google Scholar 

  55. Goldring MB, Goldring SR. Osteoarthritis. J Cell Physiol. 2007;213:626–34.

    Article  PubMed  Google Scholar 

  56. Nakamoto H, et al. Involvement of transient receptor potential vanilloid channel 2 in the induction of lubricin and suppression of ectopic endochondral ossification in mouse articular cartilage. Arthritis Rheum. 2021;73:1441–50.

    Article  Google Scholar 

  57. Nasi S, et al. The protective role of the 3-mercaptopyruvate sulfurtransferase (3-MST)-hydrogen sulfide (H2S) pathway against experimental osteoarthritis. Arthritis Res Ther. 2020;22:1–2.

    Article  Google Scholar 

  58. Grant MP, et al. Human cartilaginous endplate degeneration is induced by calcium and the extracellular calcium-sensing receptor in the intervertebral disc. Eur Cell Mater. 2016;32:137–51.

    Article  PubMed  Google Scholar 

  59. Otsuki S, et al. Extracellular sulfatases support cartilage homeostasis by regulating BMP and FGF signaling pathways. Proc Natl Acad Sci U S A. 2010;107:10202–7.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Settembre C, et al. Proteoglycan desulfation determines the efficiency of chondrocyte autophagy and the extent of FGF signaling during endochondral ossification. Genes Dev. 2008;22:2645–50.

    Article  PubMed  PubMed Central  Google Scholar 

  61. Koziel L, Kunath M, Kelly OG, Vortkamp A. Ext1-dependent heparan sulfate regulates the range of Ihh signaling during endochondral ossification. Dev Cell. 2004;6:801–13.

    Article  PubMed  Google Scholar 

  62. Sarrazin S, Lamanna WC, Esko JD. Heparan sulfate proteoglycans. Cold Spring Harb Perspect Biol. 2011;3:1–33.

    Article  Google Scholar 

  63. Guilak F, Hayes AJ, Melrose J. Perlecan in pericellular mechanosensory cell-matrix communication, extracellular matrix stabilisation and mechanoregulation of load-bearing connective tissues. Int J Mol Sci. 2021;22:1–20 Preprint at. https://doi.org/10.3390/ijms22052716.

    Article  Google Scholar 

  64. Otsuki S, Alvarez-Garcia O, Lotz MK, Neo M. Role of heparan sulfate 6-0 endosulfatases in intervertebral disc homeostasis. Histol Histopathol. 2019;34:1051–60.

    PubMed  PubMed Central  Google Scholar 

  65. Malaver-Ortega LF, Sumer H, Liu J, Verma PJ. The state of the art for pluripotent stem cells derivation in domestic ungulates. Theriogenology. 2012;78:1749–62.

    Article  PubMed  Google Scholar 

  66. Rutges JPHJ, et al. Hypertrophic differentiation and calcification during intervertebral disc degeneration. Osteoarthr Cartil. 2010;18:1487–95.

    Article  Google Scholar 

  67. Illien-Junger S, et al. AGEs induce ectopic endochondral ossification in intervertebral discs. Eur Cell Mater. 2016;32:257–70.

    Article  PubMed  PubMed Central  Google Scholar 

  68. Hallett SA, Ono W, Ono N. The hypertrophic chondrocyte: to be or not to be. Histol Histopathol. 2021;36:1021–36 Preprint at. 10.14670/HH-18-355

    PubMed  PubMed Central  Google Scholar 

  69. Rodrigues M, Kosaric N, Bonham CA, Gurtner GC. Wound healing: a cellular perspective. Physiol Rev. 2019;99:665–706.

    Article  PubMed  Google Scholar 

  70. Xie Y, et al. LMO7 is a negative feedback regulator of transforming growth factor β signaling and fibrosis. Circulation. 2019;139:679–93.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Kolle G, Georgas K, Holmes GP, Little MH, Yamada T. CRIM1, a novel gene encoding a cysteine-rich repeat protein, is developmentally regulated and implicated in vertebrate CNS development and organogenesis. Mech Dev. 2000;90:181–93.

    Article  PubMed  Google Scholar 

  72. Tang, R. et al. Gene therapy for follistatin mitigates systemic metabolic inflammation and post-traumatic arthritis in high-fat diet-induced obesity. https://www.science.org (2020).

  73. Del Monte P, et al. Effects of α-interferon on insulin-like growth factor-I, insulin-like growth factor-II and insulin-like growth factor binding protein-3 secretion by a human lung cancer cell line in vitro. J Endocrinol Investig. 2005;28:432–9.

    Article  Google Scholar 

  74. Camozzi M, et al. Identification of an antiangiogenic FGF2-binding site in the N terminus of the soluble pattern recognition receptor PTX3. J Biol Chem. 2006;281:22605–13.

    Article  PubMed  Google Scholar 

  75. Kawashima K, et al. Heparan sulfate deficiency leads to hypertrophic chondrocytes by increasing bone morphogenetic protein signaling. Osteoarthr Cartil. 2020;28:1459–70.

    Article  Google Scholar 

  76. Mundy C, et al. Synovial joint formation requires local Ext1 expression and heparan sulfate production in developing mouse embryo limbs and spine. Dev Biol. 2011;351:70–81.

    Article  PubMed  Google Scholar 

  77. Yang Y, Topol L, Lee H, Wu J, et al. Development. 2003;130:1003–15 Preprint at. https://doi.org/10.1242/dev.00324.

    Article  PubMed  Google Scholar 

  78. Liu J, et al. Wnt/β-catenin signalling: function, biological mechanisms, and therapeutic opportunities. Signal Transduct Target Ther. 2022;7:3. https://doi.org/10.1038/s41392-021-00762-6.

    Article  PubMed  PubMed Central  Google Scholar 

  79. Oichi T, Otsuru S, Usami Y, Enomoto-Iwamoto M, Iwamoto M. Wnt signaling in chondroprogenitors during long bone development and growth. Bone. 2020;137:115368. https://doi.org/10.1016/j.bone.2020.115368.

    Article  PubMed  PubMed Central  Google Scholar 

  80. Fabik J, Psutkova V, Machon O. Meis2 controls skeletal formation in the hyoid region. Front Cell Dev Biol. 2022;10:951063.

    Article  PubMed  PubMed Central  Google Scholar 

  81. Fabik J, Kovacova K, Kozmik Z, Machon O. Neural crest cells require Meis2 for patterning the mandibular arch via the Sonic hedgehog pathway. Biol Open. 2020;9:bio052043.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Selleri L, et al. Requirement for Pbx1 in skeletal patterning and programming chondrocyte proliferation and differentiation. Development. 2001;128:3543–57.

    Article  PubMed  Google Scholar 

  83. Saito T, Ikeda T, Nakamura K, Chung UI, Kawaguchi H. S100A1 and S100B, transcriptional targets of SOX trio, inhibit terminal differentiation of chondrocytes. EMBO Rep. 2007;8:504–9.

    Article  PubMed  PubMed Central  Google Scholar 

  84. Wang W, Xu J, Du B, Kirsch T. Role of the progressive ankylosis gene (ank) in cartilage mineralization. Mol Cell Biol. 2005;25:312–23.

    Article  PubMed  PubMed Central  Google Scholar 

  85. Pasold J, et al. Reduced expression of Sfrp1 during chondrogenesis and in articular chondrocytes correlates with osteoarthritis in STR/ort mice. Exp Cell Res. 2013;319:649–59.

    Article  PubMed  Google Scholar 

  86. Diaz-Romero J, et al. S100A1 and S100B expression patterns identify differentiation status of human articular chondrocytes. J Cell Physiol. 2014;229:1106–17.

    Article  PubMed  Google Scholar 

  87. Jouan Y, et al. Lin28a induces SOX9 and chondrocyte reprogramming via HMGA2 and blunts cartilage loss in mice. Sci Adv. 2022;8:eabn3106.

    Article  PubMed  PubMed Central  Google Scholar 

  88. Baird DA, et al. Identification of novel loci associated with hip shape: a meta-analysis of genomewide association studies. J Bone Miner Res. 2019;34:241–51.

    Article  PubMed  Google Scholar 

  89. Klopocki E, et al. Deletion and point mutations of PTHLH cause brachydactyly type E. Am J Hum Genet. 2010;86:434–9.

    Article  PubMed  PubMed Central  Google Scholar 

  90. Hollander JM, et al. A critical bioenergetic switch is regulated by IGF2 during murine cartilage development. Commun Biol. 2022;5:1230.

    Article  PubMed  PubMed Central  Google Scholar 

  91. Diaz-Romero J, Nesic D. S100A1 and S100B: calcium sensors at the cross-roads of multiple chondrogenic pathways. J Cell Physiol. 2017;232:1979–87 Preprint at. https://doi.org/10.1002/jcp.25720.

    Article  PubMed  Google Scholar 

  92. Chen L, Fink T, Zhang XY, Ebbesen P, Zachar V. Quantitative transcriptional profiling of ATDC5 mouse progenitor cells during chondrogenesis. Differentiation. 2005;73:350–63.

    Article  PubMed  Google Scholar 

  93. Schuurman W, et al. Zonal chondrocyte subpopulations reacquire zone-specific characteristics during in vitro redifferentiation. Am J Sports Med. 2009;37:97S–104S.

    Article  PubMed  Google Scholar 

  94. Facchini A, et al. Polyamine depletion inhibits NF-κB binding to DNA and interleukin-8 production in human chondrocytes stimulated by tumor necrosis factor-α. J Cell Physiol. 2005;204:956–63.

    Article  PubMed  PubMed Central  Google Scholar 

  95. Lao G, et al. Human tissue inhibitor of metalloproteinases-1 improved wound healing in diabetes through its anti-apoptotic effect. Exp Dermatol. 2019;28:528–35.

    Article  PubMed  Google Scholar 

  96. Madhu V, Guntur AR, Risbud MV. Role of autophagy in intervertebral disc and cartilage function: implications in health and disease. Matrix Biol. 2021;100–101:207–20.

    Article  PubMed  Google Scholar 

  97. Schafer K. The cell cycle: a review. Vet Pathol. 1998;35:461–78.

    Article  PubMed  Google Scholar 

  98. Cleverley K, et al. A novel knockout mouse for the small EDRK-rich factor 2 (Serf2) showing developmental and other deficits. Mamm Genome. 2021;32:94–103.

    Article  PubMed  PubMed Central  Google Scholar 

  99. Dumax-Vorzet A, Roboti P, High S. OST4 is a subunit of the mammalian oligosaccharyltransferase required for efficient N-glycosylation. J Cell Sci. 2013;126:2595–606.

    PubMed  PubMed Central  Google Scholar 

  100. Park KS, Bandeira E, Shelke GV, Lässer C, Lötvall J. Enhancement of therapeutic potential of mesenchymal stem cell-derived extracellular vesicles. Stem Cell Res Ther. 2019;10:1–5 Preprint at. https://doi.org/10.1186/s13287-019-1398-3.

    Article  Google Scholar 

  101. Cheng X, et al. Mesenchymal stem cells deliver exogenous miR-21 via exosomes to inhibit nucleus pulposus cell apoptosis and reduce intervertebral disc degeneration. J Cell Mol Med. 2018;22:261–76.

    Article  PubMed  Google Scholar 

  102. Piazza N, Dehghani M, Gaborski TR, Wuertz-Kozak K. Therapeutic potential of extracellular vesicles in degenerative diseases of the intervertebral disc. Front Bioeng Biotechnol. 2020;8 Preprint at https://doi.org/10.3389/fbioe.2020.00311.

  103. Kim YG, Choi J, Kim K. Mesenchymal stem cell-derived exosomes for effective cartilage tissue repair and treatment of osteoarthritis. Biotechnol J. 2020;15 Preprint at https://doi.org/10.1002/biot.202000082.

  104. Liu LT, et al. Characteristics of stem cells derived from the degenerated human intervertebral disc cartilage endplate. PLoS One. 2011;6:e26285.

    Article  PubMed  PubMed Central  Google Scholar 

  105. Luo L, et al. Cartilage endplate stem cells inhibit intervertebral disc degeneration by releasing exosomes to nucleus pulposus cells to activate Akt/autophagy. Stem Cells. 2021;39:467–81.

    Article  PubMed  Google Scholar 

  106. Luo L, et al. Cartilage endplate stem cells transdifferentiate into nucleus pulposus cells via autocrine exosomes. Front Cell Dev Biol. 2021;9:648201.

    Article  PubMed  PubMed Central  Google Scholar 

  107. Fernandes LM, et al. Single-cell RNA-seq identifies unique transcriptional landscapes of human nucleus pulposus and annulus fibrosus cells. Sci Rep. 2020;10:15263.

    Article  PubMed  PubMed Central  Google Scholar 

  108. Panebianco CJ, Dave A, Charytonowicz D, Sebra R, Iatridis JC. Single-cell RNA-sequencing atlas of bovine caudal intervertebral discs: Discovery of heterogeneous cell populations with distinct roles in homeostasis. FASEB J. 2021;35:e21919.

    Article  PubMed  Google Scholar 

  109. Wang D, et al. Single-cell transcriptomics reveals heterogeneity and intercellular crosstalk in human intervertebral disc degeneration. iScience. 2023;26:106692.

    Article  PubMed  PubMed Central  Google Scholar 

  110. Wang J, et al. Novel biomarkers of intervertebral disc cells and evidence of stem cells in the intervertebral disc. Osteoarthr Cartil. 2021;29:389–401.

    Article  Google Scholar 

  111. Cherif H, et al. Single-cell RNA-Seq analysis of cells from degenerating and non-degenerating intervertebral discs from the same individual reveals new biomarkers for intervertebral disc degeneration. Int J Mol Sci. 2022;23:3993.

    Article  PubMed  PubMed Central  Google Scholar 

  112. Gan Y, et al. Spatially defined single-cell transcriptional profiling characterizes diverse chondrocyte subtypes and nucleus pulposus progenitors in human intervertebral discs. Bone Res. 2021;9:37.

    Article  PubMed  PubMed Central  Google Scholar 

  113. Xu P, et al. Fox proteins are modular competency factors for facial cartilage and tooth specification. Development (Cambridge). 2018;145:dev165498.

    Article  Google Scholar 

  114. Xu J, et al. A Shh-Foxf-Fgf18-Shh molecular circuit regulating palate development. PLoS Genet. 2016;12:e1005769.

    Article  PubMed  PubMed Central  Google Scholar 

  115. Balic A, Adams D, Mina M. Prx1 and Prx2 cooperatively regulate the morphogenesis of the medial region of the mandibular process. Dev Dyn. 2009;238:2599–613.

    Article  PubMed  PubMed Central  Google Scholar 

  116. Yammani RR. S100 proteins in cartilage: role in arthritis. Biochim Biophys Acta Mol basis Dis. 2012;1822:600–6 Preprint at. https://doi.org/10.1016/j.bbadis.2012.01.006.

    Article  Google Scholar 

  117. Bedford L, et al. Id4 is required for the correct timing of neural differentiation. Dev Biol. 2005;280:386–95.

    Article  PubMed  Google Scholar 

  118. Dehne T, et al. Gene expression profiling of primary human articular chondrocytes in high-density micromasses reveals patterns of recovery, maintenance, re- and dedifferentiation. Gene. 2010;462:8–17.

    Article  PubMed  Google Scholar 

  119. Maeda-Uematsu A, et al. CCN2 as a novel molecule supporting energy metabolism of chondrocytes. J Cell Biochem. 2014;115:854–65.

    Article  PubMed  PubMed Central  Google Scholar 

  120. Lane RS, et al. Mitochondrial respiration and redox coupling in articular chondrocytes. Arthritis Res Ther. 2015;17:1–4.

    Article  Google Scholar 

  121. Agrawal A, et al. Normoxic stabilization of HIF-1α drives glycolytic metabolism and regulates aggrecan gene expression in nucleus pulposus cells of the rat intervertebral disk. Am J Phys Cell Phys. 2007;293:C621–31.

    Google Scholar 

  122. Kruithof PD, et al. Unraveling the role of thiosulfate sulfurtransferase in metabolic diseases. Biochim Biophys Acta Mol basis Dis. 2020;1866:165716 Preprint at. https://doi.org/10.1016/j.bbadis.2020.165716.

    Article  PubMed  Google Scholar 

  123. Ohashi Y, et al. Metabolic reprogramming in chondrocytes to promote mitochondrial respiration reduces downstream features of osteoarthritis. Sci Rep. 2021;11:15131.

    Article  PubMed  PubMed Central  Google Scholar 

  124. Hartman R, et al. Age-dependent changes in intervertebral disc cell mitochondria and bioenergetics. Eur Cell Mater. 2018;36:171–83.

    Article  PubMed  PubMed Central  Google Scholar 

  125. Mang T, et al. BMPR1A is necessary for chondrogenesis and osteogenesis, whereas BMPR1B prevents hypertrophic differentiation. J Cell Sci. 2020;133:jcs246934.

    Article  PubMed  Google Scholar 

  126. Hu J, et al. Formin 1 and filamin B physically interact to coordinate chondrocyte proliferation and differentiation in the growth plate. Hum Mol Genet. 2014;23:4663–73.

    Article  PubMed  PubMed Central  Google Scholar 

  127. Watanabe-Takano H, et al. Mechanical load regulates bone growth via periosteal Osteocrin. Cell Rep. 2021;36:109380.

    Article  PubMed  Google Scholar 

  128. Shimoda M, et al. Hyaluronan-binding protein involved in hyaluronan depolymerization controls endochondral ossification through hyaluronan metabolism. Am J Pathol. 2017;187:1162–76.

    Article  PubMed  Google Scholar 

  129. Zhu S, et al. Subchondral bone osteoclasts induce sensory innervation and osteoarthritis pain. J Clin Investig. 2019;129:1076–93.

    Article  PubMed  PubMed Central  Google Scholar 

  130. Ni S, et al. Sensory innervation in porous endplates by Netrin-1 from osteoclasts mediates PGE2-induced spinal hypersensitivity in mice. Nat Commun. 2019;10:5643.

    Article  PubMed  PubMed Central  Google Scholar 

  131. Akoum J, et al. Netrin-1 secreted by human osteoarthritic articular chondrocytes promotes angiogenesis in vitro. Cartilage. 2022;13:94–104.

    Article  PubMed  PubMed Central  Google Scholar 

  132. Loeser RF. Integrins and chondrocyte-matrix interactions in articular cartilage. Matrix Biol. 2014;39:11–6 Preprint at. https://doi.org/10.1016/j.matbio.2014.08.007.

    Article  PubMed  PubMed Central  Google Scholar 

  133. Fearing BV, Hernandez PA, Setton LA, Chahine NO. Mechanotransduction and cell biomechanics of the intervertebral disc. JOR Spine. 2018;1 Preprint at https://doi.org/10.1002/jsp2.1026.

  134. Wipff PJ, Hinz B. Integrins and the activation of latent transforming growth factor β1 - an intimate relationship. Eur J Cell Biol. 2008;87:601–15 Preprint at. https://doi.org/10.1016/j.ejcb.2008.01.012.

    Article  PubMed  Google Scholar 

  135. Zhen G, et al. Inhibition of TGF-β signaling in mesenchymal stem cells of subchondral bone attenuates osteoarthritis. Nat Med. 2013;19:704–12.

    Article  PubMed  PubMed Central  Google Scholar 

  136. Bian Q, et al. Mechanosignaling activation of TGFβ maintains intervertebral disc homeostasis. Bone Res. 2017;5:1–4.

    Article  Google Scholar 

  137. Tsutsui K, et al. ADAMTSL-6 is a novel extracellular matrix protein that binds to fibrillin-1 and promotes fibrillin-1 fibril formation. J Biol Chem. 2010;285:4870–82.

    Article  PubMed  Google Scholar 

  138. Chaudhry SS, et al. Fibrillin-1 regulates the bioavailability of TGFβ1. J Cell Biol. 2007;176:355–67.

    Article  PubMed  PubMed Central  Google Scholar 

  139. Wu Z, et al. CEP290 is essential for the initiation of ciliary transition zone assembly. PLoS Biol. 2020;18:e3001034.

    Article  PubMed  PubMed Central  Google Scholar 

  140. Dagoneau N, et al. DYNC2H1 mutations cause asphyxiating thoracic dystrophy and short rib-polydactyly syndrome, type III. Am J Hum Genet. 2009;84:706–11.

    Article  PubMed  PubMed Central  Google Scholar 

  141. Thompson CL, McFie M, Paul Chapple J, Beales P, Knight MM. Polycystin-2 is required for chondrocyte mechanotransduction and traffics to the primary cilium in response to mechanical stimulation. Int J Mol Sci. 2021;22:4313.

    Article  PubMed  PubMed Central  Google Scholar 

  142. Scott CC, Vacca F, Gruenberg J. Endosome maturation, transport and functions. Semin Cell Dev Biol. 2014;31:2–10 Preprint at. https://doi.org/10.1016/j.semcdb.2014.03.034.

    Article  PubMed  Google Scholar 

  143. Tooze SA, Abada A, Elazar Z. Endocytosis and autophagy: exploitation or cooperation? Cold Spring Harb Perspect Biol. 2014;6:a018358.

    Article  PubMed  PubMed Central  Google Scholar 

  144. Zuo R, et al. Rapamycin induced autophagy inhibits inflammation-mediated endplate degeneration by enhancing Nrf2/Keap1 signaling of cartilage endplate stem cells. Stem Cells. 2019;37:828–40.

    Article  PubMed  Google Scholar 

  145. Xu H. Autophagy protects endplate chondrocytes from intermittent cyclic mechanical tension induced calcification. Bone. 2015;75:242–3 Preprint at. https://doi.org/10.1016/j.bone.2014.10.023.

    Article  PubMed  Google Scholar 

  146. Chen K, et al. Autophagy is a protective response to the oxidative damage to endplate chondrocytes in intervertebral disc: implications for the treatment of degenerative lumbar disc. Oxidative Med Cell Longev. 2017;2017:1–9.

    Google Scholar 

  147. Jayasuriya CT, Chen Q. Potential benefits and limitations of utilizing chondroprogenitors in cell-based cartilage therapy. Connect Tissue Res. 2015;56:265–71.

    Article  PubMed  PubMed Central  Google Scholar 

  148. Choi H, Johnson ZI, Risbud MV. Understanding nucleus pulposus cell phenotype: a prerequisite for stem cell based therapies to treat intervertebral disc degeneration. Curr Stem Cell Res Ther. 2015;10:307–16.

    Article  PubMed  PubMed Central  Google Scholar 

  149. O’Brien K, Breyne K, Ughetto S, Laurent LC, Breakefield XO. RNA delivery by extracellular vesicles in mammalian cells and its applications. Nat Rev Mol Cell Biol. 2020;21:585–606 Preprint at. https://doi.org/10.1038/s41580-020-0251-y.

    Article  PubMed  PubMed Central  Google Scholar 

  150. Zhang Y, Liu Y, Liu H, Tang WH. Exosomes: biogenesis, biologic function and clinical potential. Cell Biosci. 2019;9:1–8.

    Article  Google Scholar 

  151. Van Den Akker GGH, Caron MMJ, Peffers MJ, Welting TJM. Ribosome dysfunction in osteoarthritis. Curr Opin Rheumatol. 2022;34:61–7 Preprint at. https://doi.org/10.1097/BOR.0000000000000858.

    Article  PubMed  Google Scholar 

  152. Peffers MJ, et al. SnoRNA signatures in cartilage ageing and osteoarthritis. Sci Rep. 2020;10:10641.

    Article  PubMed  PubMed Central  Google Scholar 

  153. Pimlott Z, et al. Small nucleolar RNAs as mediators of oxidative stress in cross species cartilage and osteoarthritis. Osteoarthr Cartil. 2020;28:S342.

    Article  Google Scholar 

  154. Chabronova A, et al. Ribosomal RNA-based epitranscriptomic regulation of chondrocyte translation and proteome in osteoarthritis. Osteoarthr Cartil. 2023;31:374–85.

    Article  Google Scholar 

  155. van den Akker GG, et al. TNF-alpha induces FGF1 IRES mediated messenger RNA translation in chondrocytes. Osteoarthr Cartil. 2021;29:S50–1.

    Article  Google Scholar 

  156. Komar AA, Hatzoglou M. Cellular IRES-mediated translation: The war of ITAFs in pathophysiological states. Cell Cycle. 2011;10:229–40 Preprint at. https://doi.org/10.4161/cc.10.2.14472.

    Article  PubMed  PubMed Central  Google Scholar 

  157. Marques R, Lacerda R, Romão L. Internal ribosome entry site (IRES)-mediated translation and its potential for novel mRNA-based therapy development. Biomedicines. 2022;10:1865 Preprint at. https://doi.org/10.3390/biomedicines10081865.

    Article  PubMed  PubMed Central  Google Scholar 

  158. Yang Z, et al. Dysregulated COL3A1 and RPL8, RPS16, and RPS23 in disc degeneration revealed by bioinformatics methods. Spine (Phila Pa 1976). 2015;40:E745–51.

    Article  PubMed  Google Scholar 

  159. Holm S, Selstam G. Oxygen tension alterations in the intervertebral disc as a response to changes in the arterial blood. Ups J Med Sci. 1982;87:163–74.

    Article  PubMed  Google Scholar 

  160. Amorim JA, et al. Mitochondrial and metabolic dysfunction in ageing and age-related diseases. Nat Rev Endocrinol. 2022;18:243–58.

    Article  PubMed  PubMed Central  Google Scholar 

  161. Saberi M, Zhang X, Mobasheri A. Targeting mitochondrial dysfunction with small molecules in intervertebral disc aging and degeneration. Geroscience. 2021;43:517–37.

    Article  PubMed  PubMed Central  Google Scholar 

  162. Cisewski SE, et al. Comparison of oxygen consumption rates of nondegenerate and degenerate human intervertebral disc cells. Spine (Phila Pa 1976). 2018;43:E60–7.

    Article  PubMed  Google Scholar 

  163. Ma B, et al. Gene expression profiling of dedifferentiated human articular chondrocytes in monolayer culture. Osteoarthr Cartil. 2013;21:599–603.

    Article  Google Scholar 

  164. Le Maitre CL, et al. Altered integrin mechanotransduction in human nucleus pulposus cells derived from degenerated discs. Arthritis Rheum. 2009;60:460–9.

    Article  PubMed  Google Scholar 

  165. Benneker LM, Heini PF, Anderson SE, Alini M, Ito K. Correlation of radiographic and MRI parameters to morphological and biochemical assessment of intervertebral disc degeneration. Eur Spine J. 2005;14:27–35.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

Bulk- and single-cell RNA sequencing was performed by The Ohio State University Comprehensive Cancer Center’s Genomics Shared Resource supported by NIH grant P30-CA016058.

Funding

Funding for this work was supported by the award UL1TR002733 from the National Center for Advancing Translational Sciences and NSF CAREER grant 2143123.

Author information

Authors and Affiliations

Authors

Contributions

K.K., S.T., S.K., B.W., and D.P. designed the experiments. K.K. and S.T. collected and prepared cells for RNA-Seq. P.S., A.H., and M.P. performed the bioinformatic analysis. K.K., S.T., W.X., and K.D. performed the literature review. K.K., P.S., A.H., and M.P. prepared the figures and tables. K.K., S.T., P.S., M.P., B.W., and D.P. contributed to the data interpretation and analysis. K.K. drafted the manuscript with feedback from co-authors. All authors have read and approved the final manuscript draft for submission.

Corresponding author

Correspondence to Devina Purmessur.

Ethics declarations

Ethics approval and consent to participate

All cadaveric human tissue samples were obtained through Ohio State University’s Cooperative Human Tissue Network and were obtained through an Institutional Review Board exemption as this research is classified as non-human subject research.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

Number of publications, search terms, and PubMed IDs for all identified papers for all genes in each tissue of interest.

Additional file 2: Table S2.

Fold change and adjusted p-values for all significantly D.E. genes. Table S3. General breakdown of genes with prior publications in cartilaginous cells (N=31 total genes, q<0.05). Table S4. General breakdown of genes with no prior publications in cartilaginous cells (N=45 total genes, q<0.05). Table S5. Screening of cell types and corresponding markers possibly isolated in CEP samples. Table S6. Cell types and corresponding markers identified in current single-cell RNA-Sequencing analysis.

Additional file 3: Figure S1.

Pre-filtering quality control. Violin plots showing RNA unique features, total RNA count, mitochondrial RNA content, ribosomal RNA content, hemoglobin RNA content, and platelet RNA content for all cells are shown.

Additional file 4: Table S7.

Percentage expressed and scaled expressed of differentially expressed genes from bulk RNA sequencing analysis, by cluster. 

Additional file 5: Table S8.

Functional enrichment values and gene lists for Gene Ontology pathways, in chondrocytes.

Additional file 6: Table S9.

Functional enrichment values and gene lists for Gene Ontology pathways, in progenitor cells.

Additional file 7: Figure S2.

Dot plot of top 3 markers of each cluster identified in single-cell RNA-Sequencing analysis.

Additional file 8: Figure S3.

Dot plots of markers associated with different phases of mitosis. Markers were selected from the following gene sets from https://maayanlab.cloud/Harmonizome/: 1) Prophase (GeneRIF Biological Term Annotations), 2) Prometaphase (GeneRIF Biological Term Annotations), and 3) Metaphase (GeneRIF Biological Term Annotations).

Additional file 9: Figure S4.

Dot plots of markers associated with different phases of mitosis. Markers were selected from the following genes sets from https://maayanlab.cloud/Harmonizome/: 1) Anaphase (GeneRIF Biological Term Annotations), 2) Telophase (GeneRIF Biological Term Annotations), and 3) Cytokinesis (GO Biological Process Annotations).

Additional file 10: Figure S5.

Module score uMAP of mitosis-related genes from Fig. S3 and Fig. S4.

Additional file 11: Figure S6.

Dot plot of markers associated with G1/S phase of cell cycle. Markers were selected from gene sets from https://maayanlab.cloud/Harmonizome/.

Additional file 12: Figure S7.

Module score uMAP of G1/S transition genes from Fig. S6.

Additional file 13: Figure S8.

Gene ontology comparing Non-degenerated and Degenerated Chondrocyte 1. Red boxes indicate pathways associated with ribosomes, protein translation, and mitochondrial function. “Activated” pathways are enriched in the degenerated sample and “Suppressed” pathways are enriched in the non-degenerated sample.

Additional file 14: Figure S9.

Gene ontology comparing Non-degenerated and Degenerated Chondrocyte 2. Red boxes indicate pathways associated with ribosomes, protein translation, and mitochondrial function. “Activated” pathways are enriched in the degenerated sample and “Suppressed” pathways are enriched in the non-degenerated sample.

Additional file 15: Figure S10.

Gene ontology comparing Non-degenerated and Degenerated Chondrocyte 3. Red boxes indicate pathways associated with ribosomes, protein translation, and mitochondrial function. “Activated” pathways are enriched in the degenerated sample and “Suppressed” pathways are enriched in the non-degenerated sample.

Additional file 16: Figure S11.

Gene ontology comparing Non-degenerated and Degenerated Chondroprogenitors. Red boxes indicate pathways associated with ribosomes, protein translation, and mitochondrial function. “Activated” pathways are enriched in the degenerated sample and “Suppressed” pathways are enriched in the non-degenerated sample.

Additional file 17: Figure S12.

Gene ontology comparing Non-degenerated and Degenerated MSCs. Red boxes indicate pathways associated with ribosomes, protein translation, and mitochondrial function. “Activated” pathways are enriched in the degenerated sample and “Suppressed” pathways are enriched in the non-degenerated sample.

Additional file 18: Figure S13.

Gene ontology comparing Non-degenerated and Degenerated Proliferating MSCs. Red boxes indicate pathways associated with ribosomes, protein translation, and mitochondrial function. “Activated” pathways are enriched in the degenerated sample and “Suppressed” pathways are enriched in the non-degenerated sample.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kuchynsky, K., Stevens, P., Hite, A. et al. Transcriptional profiling of human cartilage endplate cells identifies novel genes and cell clusters underlying degenerated and non-degenerated phenotypes. Arthritis Res Ther 26, 12 (2024). https://doi.org/10.1186/s13075-023-03220-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13075-023-03220-6

Keywords