Single-cell RNA-seq identifies unique transcriptional landscapes of human nucleus pulposus and annulus fibrosus cells

Abstract

Intervertebral disc (IVD) disease (IDD) is a complex, multifactorial disease. While various aspects of IDD progression have been reported, the underlying molecular pathways and transcriptional networks that govern the maintenance of healthy nucleus pulposus (NP) and annulus fibrosus (AF) have not been fully elucidated. We defined the transcriptome map of healthy human IVD by performing single-cell RNA-sequencing (scRNA-seq) in primary AF and NP cells isolated from non-degenerated lumbar disc. Our systematic and comprehensive analyses revealed distinct genetic architecture of human NP and AF compartments and identified 2,196 differentially expressed genes. Gene enrichment analysis showed that SFRP1, BIRC5, CYTL1, ESM1 and CCNB2 genes were highly expressed in the AF cells; whereas, COL2A1, DSC3, COL9A3, COL11A1, and ANGPTL7 were mostly expressed in the NP cells. Further, functional annotation clustering analysis revealed the enrichment of receptor signaling pathways genes in AF cells, while NP cells showed high expression of genes related to the protein synthesis machinery. Subsequent interaction network analysis revealed a structured network of extracellular matrix genes in NP compartments. Our regulatory network analysis identified FOXM1 and KDM4E as signature transcription factor of AF and NP respectively, which might be involved in the regulation of core genes of AF and NP transcriptome.

Figure 1

Enrichment of tissue-specific markers and distinct segregation of genes in NP and AF cells. (a) Schematic representation of SCT workflow from tissue acquisition to analysis. The IVD image was digitally captured using KodakEasyShare v550 camera and inner dashed circle were traced using MS-Power Point drawing tool to demarcate NP and AF. (b) Relative expression of COL2A1 shows enrichment of the NP marker gene in NP Cells compared to AF cells. (c) Relative Expression of CYTL1 shows higher expression of the AF marker gene in AF cells compared to NP cells. (d) tSNE plot showing clear segregation of primary NP and AF cells. The tSNE plot was genereated using Seurat package in R version 3.0.

Figure 2

Differential Expression of genes in AF and NP cells. (a) Volcano plot depicting differentially expressed genes in AF and NP cells. Red dots represent genes expressed at higher levels in AF cells while blue dots represent genes with higher expression levelsin NP cells. Y-axis denotes − log10 P values while X-axis shows log2 fold change values. Volcano plot was generated using GraphPad Prism version 8.2.0. (b) Relative gene expression of LRRC17, AK5, SFRP1 and KIAA0101 in AF cells from 3 human samples expressed as a percentage of expression in NP. (c) Relative gene expression of COL11A1, DSC3, COL9A3 and FAM46B genes in NP cells from 3 human samples expressed as a percentage of expression in AF cells. Black, Pink and Blue bars represent gene expression of cells from 24-year-old female, 35-year-old male and 18-year-old female respectively. Bar diagram for gene expression was generated using GraphPad Prism version 8.2.0.

Figure 3

Single cell trancriptomics (SCT) analysis identifies novel genes in AF and NP cells: (a) Venn diagram depicting common and unique genes detected by microarray and SCT analysis. (b) Violin plot showing genome-wide similarity between microarray and SCT analysis data set. (c) Scatter plot displaying the similarity in the expression profile of the 214 common genes between SCT and Microarray data. Each dot represents expression value of a gene in the data set. All the plots were generated using GraphPad Prism version 8.2.0.

Figure 4

Gene Ontology analysis reveal compartment-specific enrichment of GO pathways in AF and NP cells. (a) Advanced bubble plot depicting molecular function-related pathways enriched in AF cells. (b) Advanced bubble plot showing molecular function-related pathways enriched in NP cells. Y-axis Labels represent pathway names while and enrichment score is shown on X-axis. The size and color of the bubble represent the number of genes enriched in each pathway and their enrichment significance. The bubble plots were generated using GraphPad Prism version 8.2.0. (c) Metabolic pathway analysis of glycolysis, TCA cycle, electron transport chain and pyruvate metabolism pathways in AF. (d) Metabolic Pathway analysis of glycolysis, TCA cycle, electron transport chain and pyruvate metabolism pathways in NP cells. Heatmap shows log10 transformed fold change in gene expression and were generated using GraphPad Prism version 8.2.0.

Figure 5

Core matrisome analysis of AF and NP cells show differential expression of genes. (a) Venn diagram depicting the total Core matrisome genes, the Core matrisome genes that are present in the IVD, and the total number of differentially expressed genes in our analysis. (b) Venn diagram representing the total Core matrisome associated genes, the Core matrisome associated genes that we detect in the IVD and the total number of differentially expressed genes in our analysis. (c) Pie chart showing the distribution of Core matrisome genes in AF and NP (AF is shown in red and NP in blue). (d) Pie chart showing the distribution of Core matrisome associated genes in AF and NP (AF is shown in red and NP in blue). (e) Heatmap showing the differential gene expression profile of Core matrisome genes in human NP cells. (f) Heatmap showing the differential gene expression profile of Core matrisome associated genes in human NP cells. (g) Heatmap representing the differential gene expression profile of Core matrisome genes in human AF cells. (h) Heatmap representing the differential gene expression profile of Core matrisome associated genes in human AF cells. All heatmaps were generated using GraphPad Prism version 8.2.0.

Figure 6

Regulatory network analysis identifies AF and NP-specific TF-targeted regulatory networks. (a) AF specific network cluster with 26 genes and 311 interactions. Red lines represent interactions, and yellow circles represent genes. (b) NP-specific network cluster with 7 Genes and 20 interactions. Red lines represent interactions, and yellow circles represent genes. The network clusters were made using Cytoscape (version: 3.7.1). (c) TF-targeted regulatory network showing FOXM1 and the genes regulated by it in the AF. (d) TF-targeted regulatory network showing KDM4E and the genes regulated by it in the NP. The regulatory networks were generated by iRegulon (version 1.3) package using Cytoscape (version: 3.7.1). (e) Relative expression of FOXM1 expressed as a percentage of expression in NP cells. (f) Relative expression of KDM4E expressed as a percentage of expression in AF cells. Bar diagram for gene expression was generated using GraphPad Prism version 8.2.0.