This article is mentioned in:

Introduction

Intervertebral disc degeneration (IDD) is characterized by damage of the disc structure, progressive loss of water content and proteoglycan in the extracellular matrix (ECM) (1). IDD is considered a predominant source of spine-related diseases and chronic lower back pain, which causes a major economic and social burden affecting millions of people worldwide (2). Current treatment options for IDD, including discectomy (3), intradiscal electrothermal therapy (4) and arthroplasty (5), only address the clinical symptoms of IDD and remain limited with unpredictable outcomes (6). Therefore, understanding the pathophysiology and molecular mechanism underlying IDD appears to be imperative for diagnosis and developing novel therapeutic approaches.

The intervertebral disc consists of three morphologically distinct regions, which are the nucleus pulposus (NP), annulus fibrosis and cartilaginous endplates (7). The central NP is a gelatinous matrix that is composed of large aggregating proteoglycans and a loose network of collagen. The peripheral annulus fibrosis encases the nucleus pulposus and is rich in type I collagen. The cartilaginous endplates contain the peripheral vasculature, which can nourish the disc (1). The main morphological manifestations of IDD are vertebral instability, disc herniation and spinal stenosis (8). Significant changes in morphology, structure and composition are accompanied by specific changes in the disc with aging and degeneration, including alteration of the elastic modulus and swelling pressure of the nucleus pulposus (9). Recently, histological grading schemes for assessing human IDD have been developed (10,11). A 5-level grading system for lumbar disc degeneration, proposed by Pfirrmann et al (11), was developed according to T2-weighted magnetic resonance images. Additionally, the Thompson scoring system scores disc degeneration over the spectrum from healthier discs (grades I and II) to advanced degeneration (grade V, the most advanced stage of degeneration) (10). However, there is a lack of studies on the associations of the grades of disc degeneration.

IDD is a complex multi-factorial process (12). Several factors, such as biological, mechanical and genetic factors, are widely considered as significant contributors to the disc degenerative process (1). For example, Bachmeier et al (13) identified that matrix metalloproteinase-3 (MMP-3), which was essential for matrix degeneration, had an essential role in lumbar disc herniation and degeneration. Additionally, Takahashi et al (14) revealed that the polymorphism 5A allele, which often occurs in the promoter region of the gene that regulates MMP-3 expression, was a possible risk factor for accelerated IDD in the elderly. In addition, Pratsinis et al (15) reported that platelet-derived growth factor, insulin-like growth factor-I and basic fibroblast growth factor could stimulate the proliferation of intervertebral disc cells via the activation of the extracellular-signal regulated kinase and Akt signaling pathways. By contrast, microRNAs (miRNAs) that can regulate RNA degradation or repression of translation have been identified as key regulators in numerous biological processes (16). Recently, Yu et al (17) demonstrated that miR-10b could promote nucleus pulposus cell proliferation by targeting homeobox D10 through the ras homolog family member C-Akt pathway in IDD. However, the molecular mechanism of IDD is not fully understood and further investigations are required.

Recently, Chen et al (18) reported that mitogen-activated protein (MAP) kinase kinase 6 (MAP2K6) and Ras homologous-related BTB domain containing 2 may have significant roles in the progression of grade III and IV disc degeneration, respectively. Additionally, Tang et al (8) demonstrated that differentially expressed genes (DEGs) identified in degenerative disc tissue samples were mainly associated with transforming growth factor-β and the ECM. Furthermore, Wang et al (19) reported that deregulated miR-155 could promote Fas-mediated apoptosis in human IDD by targeting Fas-associated death domain-containing protein (FADD) and caspase-3. In the present study, 3 microarray datasets were integrated and micro-array analysis was used to identify the DEGs, grade-specific genes and differentially expressed miRNAs in the IDD samples compared with their corresponding controls. In addition, comprehensive bioinformatics was used to analyze the significant functions and pathways, and to analyze the miRNA-target gene-pathway interaction associations. The study aimed to identify the significant genes and miRNA alternations in the progression of IDD and gain more insights into the molecular mechanisms of disc degeneration. Understanding these mechanisms can aid in exploring an appropriate molecular target and developing new therapeutic methods for IDD.

Materials and methods

Microarray data

Microarray datasets of GSE19943 (19), GSE15227 (20) and GSE34095 (21) were downloaded from the National Center of Biotechnology Information (NCBI), Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The miRNA expression profile GSE19943 included 3 degenerative NP samples that were collected from patients with IDD and 3 NP controls derived from patients with scoliosis, which was based on the platform GPL9946 Exiqon human miRCURY LNA™ microRNA Array V11.0. The GSE15227 dataset included 5 grade II, 7 grade III, and 3 grade IV discs from patients with herniated discs and IDD that were scored using the Thompson scoring system, which was based on the platform GPL1352 [U133_X3P] Affymetrix Human X3P array. The gene expression profile GSE34095 included 3 degenerative NP samples that were harvested from elderly patients with IDD and 3 non-degenerative samples derived from younger patients with adolescent idiopathic scoliosis as controls, which was based on the platform GPL96 [HG-U133A] Affymetrix Human Genome U133A array.

DEGs and differentially expressed miRNAs screening

The t-test in the limma package (22) was used to identify DEGs in degenerative samples compared with controls in GSE34095, DEGs in 7 grade III and 3 grade IV samples compared with 5 grade II samples in GSE15227, and to identify differentially expressed miRNAs in 3 degenerated samples compared with controls in GSE19943. The false discovery rate (FDR) was applied to perform multiple testing corrections using the Benjamini and Hochberg method (23). The threshold for the DEGs and the differentially expressed miRNAs were set as FDR <0.05.

miRNA-target gene and KEGG pathway data collection

Human miRNAs and their associated targets were downloaded from 3 high-quality online miRNA reference databases, which store manual collections of experimentally supported miRNA targets, miRecords (24), miRTarBase (25) and TarBase 6.0 (26). In addition, all pathways and associated genes were downloaded from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (27).

Grade-specific network construction

DEGs identified in grade III and IV samples were imported into Cytoscape software (28) to create protein-protein interaction (PPI) network visualizations, respectively. The source of the interaction network database was the Search Tool of the Retrieval of Interacting Genes 9.0 database, which is a comprehensive database containing functional links between proteins that are experimentally derived, as well as links predicted by text mining and comparative genomics (29). The threshold was set as confidence score, 0.4. Furthermore, the specific genes that were only enriched in one of the grade-specific networks and the common genes of the 2 grade-specific networks were extracted. The parameters describing the network topology were calculated using the Cytoscape plugin Network Analyzer (30). Logistic regression analysis (31) was conducted by SPSS 19.0 (IBM, Corp., Armonk, NY, USA) to select the most noteworthy network features.

Cluster analysis

To test the availability of these 2 grade-specific networks to reveal different stages of IDD progression, clustering analysis was performed on all genes in the networks and the top 20% nodes that had a higher degree in the network using Cytoscape plugin Network Analyzer (30), respectively. The results are represented by heat-maps.

GO and pathway enrichment analysis of DEGs and the specific genes

GO and KEGG enrichment analysis were performed for the specific genes using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (32). The threshold was set as P<0.05.

Furthermore, KEGG enrichment analyses were also performed for DEGs that were identified based on GSE34095 using DAVID. The threshold was set as P<0.05.

Measurement of miRNA-pathway interactions

With the identified differentially expressed miRNAs, as well as their target genes and the genes enriched in the KEGG pathways, Fisher's exact test (33) was used to assess the corresponding significance of the miRNA pathway.

Fisher's exact test was based on the hypergeometric distribution to combine the results of the proportion of the miRNA target gene in the functional gene set and the proportion of miRNA target gene in the whole genome. A P-value for the null hypothesis of the Fisher's exact test was examined as the genes that belong to the target genes of miRNAs and the genes that belong to the functional genes enriched.

The resulting P-value depicts the probability that the examined pathway is significantly enriched with gene targets of the selected miRNAs, the probability that at least x functional genes are enriched in the K target genes of miRNA. The P-value can be expressed as:

Results

Screening of DEGs and differentially expressed miRNAs

For GSE15227, 846 and 1,137 DEGs were identified in grade III and IV discs, respectively. For GSE34095, a total of 961 DEGs were identified in the IDD samples compared with the controls. Furthermore, for GSE19943, 77 differentially expressed miRNAs were identified in the degenerative NP samples compared with the controls.

Grade-specific network construction

The grade-specific networks are shown in Figs. 1 and 2. In the present study, 746 grade III-specific genes and 964 grade IV-specific genes were identified. Tumor protein p53 (TP53) was a hub gene in the grade III-specific network and ubiquitin C (UBC) was identified to be a hub gene in the grade IV-specific network. Additionally, there were 78 common genes in these 2 grade-specific networks. Network topology analysis showed that a total of 16 network features were identified. In addition, statistical regression analysis revealed that 6 significant network features were obtained, as shown in Table I, including average shortest path length, betweenness centrality, closeness centrality, neighborhood connectivity, radiality and stress.

Figure 1 Grade III-specific network.

Figure 2 Grade IV-specific network.

Table I Logistic regression.

Table I

Logistic regression.

Feature	P-value	Exp (B)
Average shortest path length	<0.0001	0.050
Betweenness centrality	<0.0001	<0.0001
Closeness centrality	<0.0001	<0.0001
Neighborhood connectivity	0.030	1.015
Radiality	0.036	<0.0001
Stress	<0.0001	1.010
Constant	0.001	1.08×1011

[i] Exp (B) represents odd ratio in the logistic regression model.

Clustering analysis

The clustering analysis was used to group the genes and samples on the basis of similarities of gene expression. The results of clustering analysis using all the specific genes in these 2 networks are shown in Fig. 3A. The result showed that 2 grade III samples (GSM380301 and GSM380305) were grouped into the region of grade IV samples; however, the correlation was lower compared to the correlation of the grade IV samples. By contrast, the results of clustering analysis using the top 20% genes with a higher degree in these 2 networks are shown in Fig. 3B. Consistently, the 2 grade III samples, GSM380301 and GSM380305, were grouped into the region of grade IV samples.

Figure 3 Heat-maps of specific genes. (A) Heat-maps of all the specific genes in 2 grade-specific networks. (B) Heat-maps of the top 20% specific genes that had higher degrees in the grade-specific networks. The red color indicates a high expression level and the green color indicates a low expression level of a specific gene in a particular specimen.

GO and pathway enrichment analysis of specific genes

GO and pathway analysis indicated that grade III- and IV-specific genes were significantly enriched in different GO terms and KEGG pathways. A total of 23 GO terms for grade III-specific genes were enriched, including angiogenesis, adult walking behavior and positive regulation of kinase activity (Table II). GO terms for grade IV-specific genes were mainly associated with the regulation of ubiquitin-protein ligase activity, such as positive regulation of ubiquitin-protein ligase activity, and regulation of protein ubiquitination (Table III). In addition, 4 KEGG pathways for grade III-specific genes were significantly enriched, including Alzheimer's disease, oxidative phosphorylation, Huntington's disease and Parkinson's disease. A total of 5 KEGG pathways for grade IV-specific genes were significantly enriched, including viral myocarditis, graft-versus-host disease, type I diabetes mellitus, allograft rejection and cell adhesion molecules (Table IV).

Table II GO terms of grade III-specific genes.

Table II

GO terms of grade III-specific genes.

GO term	P-value
Angiogenesis	0.007
Adult walking behavior	0.007
Retina development in camera-type eye	0.008
Activation of protein kinase activity	0.008
Response to DNA damage stimulus	0.008
Negative regulation of biosynthetic process	0.009
Negative regulation of cellular biosynthetic process	0.011
Negative regulation of macromolecule biosynthetic process	0.014
Positive regulation of protein kinase activity	0.015
Positive regulation of kinase activity	0.020
Cytokine secretion	0.024
Negative regulation of gene expression	0.024
Negative regulation of transcription	0.025
Positive regulation of transferase activity	0.026
Negative regulation of nucleobase, nucleoside, nucleotide and nucleic acid metabolic process	0.029
Negative regulation of nitrogen compound metabolic process	0.034
Vasculature development	0.037
Patterning of blood vessels	0.038
Tube morphogenesis	0.039
Regulation of protein kinase activity	0.043
Adult locomotory behavior	0.044
Actin filament organization	0.044
Cell-substrate junction assembly	0.048
Negative regulation of macromolecule metabolic process	0.050

[i] GO, Gene Ontology.

Table III GO terms of grade IV-specific genes.

Table III

GO terms of grade IV-specific genes.

GO terms	P-value
Cellular macromolecule localization	0.036
Cellular protein localization	0.030
Establishment of protein localization	0.026
Intracellular protein transport	0.038
Positive regulation of ligase activity	0.031
Positive regulation of protein ubiquitination	0.039
Positive regulation of ubiquitin-protein ligase activity	0.019
Positive regulation of ubiquitin-protein ligase activity during mitotic cell cycle	0.013
Protein transport	0.036
Regulation of ligase activity	0.025
Regulation of protein ubiquitination	0.022
Regulation of ubiquitin-protein ligase activity	0.015
Regulation of ubiquitin-protein ligase activity during mitotic cell cycle	0.022

[i] GO, Gene Ontology.

Table IV Pathway enrichment of grade III- and IV-specific genes.

Table IV

Pathway enrichment of grade III- and IV-specific genes.

Pathway	P-value
Grade III
hsa05010: Alzheimer's disease	0.009
hsa00190: Oxidative phosphorylation	0.011
hsa05016: Huntington's disease	0.012
hsa05012: Parkinson's disease	0.021
Grade IV
hsa05416: Viral myocarditis	9.37×10−4
hsa05332: Graft-versus-host disease	0.012
hsa04940: Type I diabetes mellitus	0.017
hsa05330: Allograft rejection	0.040
hsa04514: Cell adhesion molecules	0.049

miRNA-pathway-genes complex regulatory associations analysis

Subsequent to combining miRecords, miRTarBase and TarBase 6.0 databases, 5,489 miRNA-target pairs were collected in the study, including 482 miRNAs and 2,331 target genes. miRNA-pathway interactions analysis showed that there were 217 interaction pairs involved in 35 significant pathways. Furthermore, the pathway enrichment analysis for DEGs in the degenerative NP samples compared with the controls showed that 21 KEGG pathways were significantly enriched. There were 6 pathways also regulated by miRNAs, as shown in Table V, including focal adhesion, the ErbB signaling pathway, the calcium signaling pathway, the MAPK signaling pathway, apoptosis and pathways in cancer. Additionally, only miR-129-5p was involved in the calcium signaling pathway, apoptosis and pathways in cancer.

Discussion

In the present study, microarray analysis showed that a total of 846 and 1,137 DEGs were identified in grade III and IV discs, respectively. In addition, 961 DEGs were identified in the IDD samples in GSE34095 and 77 differentially expressed miRNAs were identified in the IDD samples. TP53 was a hub gene in the grade III-specific network and UBC was identified as a hub gene in grade IV-specific network. Six significant features were identified by grade-specific network topology analysis. Differentially expressed miRNAs were identified to participate in 35 pathways, among which 6 pathways were significantly enriched by DEGs, including apoptosis.

The general concepts, such as centrality, communicability and betweenness, quantify the important features in a network (34). Estrada (35) demonstrated that subgraph centrality could be applied to the identification of essential proteins in PPI networks. Additionally, the threshold for inclusion was P<0.05 (significant) in the logistic regression model (36). In line with the previous study, the present study showed that 6 network features (average shortest path length, betweenness centrality, closeness centrality, neighborhood connectivity, radiality and stress) had significant differences in grade III- and IV-specific networks (Table I). On this basis, these 6 features may have a significant association with the grade of discs degeneration and could be useful for assessing the IDD grade.

TP53 is a tumor-suppressor protein containing DNA binding, transcriptional activation and oligomerization domains (37). Vaghefi et al (38) revealed that deacetylation of p53 was a nerve growth factor (NGF)-dependent post-translational mechanism of p53 activation. Additionally, Richardson et al (39) reported that NGF had increased expression in the painful degenerate intervertebral disc. In addition, Liu et al (40) demonstrated that VEGF and p53 were working simultaneously in neovascularization and infiltration, and accelerating rat IDD progression. In line with the previous study, TP53 was identified as a hub gene in the grade III-specific network, which was derived from the microarray data of GSE15227. Furthermore, specimens of GSE15227 are annulus fibrosus cells in cultures from patients with herniated discs and IDD. Thus, we hypothesized that TP53 may participate in neovascularization and infiltration in annulus fibrosus in IDD progression, which could also help the IDD grading system.

Additionally, the present study identified that the GO terms for grade IV-specific genes were mainly associated with the regulation of ubiquitin-protein ligase activity and UBC was a hub gene in the grade IV-specific network derived from the data of GSE15227. UBC encodes a poly-ubiquitin precursor that can conjugate with different residues and lead to various effects within a cell (41). Yew (42) demonstrated that G1- and S-phase events in vertebrates, which were essential for cell proliferation, were specifically mediated by a multitude of regulators associated with ubiquitination, of which UBC was required. Furthermore, Gruber et al (43) reported that cells in IDD lost their ability to proliferate and were subject to senescence. Thus, we hypothesized that UBC may have a crucial role in inhibiting cell proliferation of annulus fibrosus in IDD progression and could be used for grading disc degeneration.

Furthermore, the present study showed that miR-129-5p could participate in different pathways, including apoptosis. The study by Kohyama et al (44) reported that induction of the apoptosis of disc cells had an important role in the pathogenesis of disc degeneration. In addition, Li et al (45) revealed that downregulation of miR-129-5p could inhibit cell growth and induce apoptosis in laryngeal squamous cell carcinoma by targeting adenomatous polyposis coli. In the present study, the miRNA-pathway interaction analysis was based on the information of differentially expressed miRNAs derived from GSE19943 and the pathway enrichment analysis for DEGs in degenerative NP samples was derived from GSE34095. Additionally, tissues of GSE19943 and GSE34095 are NP from patients with IDD and scoliosis. Furthermore, 6 overlapping pathways were identified according to the results of miRNA-pathway interactions analysis and pathway enrichment analysis for DEGs, including apoptosis. Therefore, we hypothesized that miR-129-5p may have a crucial role in modulating cell apoptosis in NP in IDD progression and miR-129-5p could be a potential therapeutic target for IDD. However, more experiments and further investigations are required to verify this finding.

In conclusion, the present study identified that key genes (TP53 and UBC) and miR-129-5p may participate in the mechanism of IDD progression. Significant network features identified in the study may aid in assessing the grade of disc degeneration. TP53 may have an essential role in neovascularization and infiltration in annulus fibrosus in grade III of IDD progression and could help the IDD grading system. Additionally, UBC may have a crucial role in inhibiting cell proliferation of annulus fibrosus in IDD progression and could be used to grade disc degeneration. Furthermore, miR-129-5p may have a crucial role in modulating cell apoptosis in NP in IDD progression and miR-129-5p could be a potential therapeutic target for IDD. However, further experiments and studies are required to confirm these results.

Acknowledgments

The present study was supported by the National Natural Science Foundation of China: The effect and mechanism of cytoskeletal elements in the mechanotransduction pathway within intervertebral disc cells and its role in intervertebral disc degeneration (grant no. 81471131) and the National Natural Science Foundation of China: The effect and mechanism of chondroitin sulphate proteoglycans in stem cell niche during stem cell differentiation (grant no. 31300675).

References