Identification of Potential Biomarkers for Intervertebral Disc Degeneration Using the Genome-Wide Expression Analysis

Abstract

Intervertebral disc degeneration (IDD) is the major cause of low back pain. The current study was aimed to further elucidate the mechanisms underlying it. Microarray data sets GSE70362 containing Thompson degeneration grades I–V were divided into the control and the degenerative group and were analyzed. Differentially expressed genes (DEGs) were screened and clustered, followed by functional enrichment analysis. Then the protein–protein interaction (PPI) network and the microRNA (miRNA) regulatory network were constructed and integratedly analyzed. Finally, modules from the integrated network were mined, and gene ontology and pathway enrichment analysis were performed. DEGs in three clusters had the overall expression trend with the Thompson grades. The upregulated DEGs were associated with protein transport and localization, while the downregulated DEGs were enriched in membrane lipid metabolic process and endocytosis. After the integrated analysis of PPI and miRNA–target interactions, some hub genes such as HSP90B1, RPS4Y1, RPL15, and UTY, as well as hub miRNAs, including miR-124a and miR-506, were screened. Finally, modules in the integrated network were functionally associated with protein targeting, peptide processing, and metabolic process, as well as protein folding. Taken together, our data showed that the protein synthesis, targeting, and localization can be greatly changed to contribute to the progress of IDD. Besides, the related genes such as RPS4Y1 and HSP90B1 would be used as diagnostic and therapeutic targets for IDD.

Highlights

Microarray data sets GSE70362 containing degeneration grades I–V were analyzed.

Protein synthesis, targeting, and localization may contribute to degenerated intervertebral disc.

RPS4Y1 and HSP90B1 would be used as diagnostic and therapeutic targets for the disease.

1. Introduction

Intervertebral disc degeneration (IDD), one of the causes of lower back pain, has led to a significant economic and social burden worldwide (Raj, 2008; Millecamps et al., 2012; Jensen et al., 2017). It is well known that the IDD occurs with age, and it could be caused by a variety of factors, such as damage, genetic, apoptosis, microenvironmental changes, and inflammation (Zhang et al., 2016). Moreover, IDD has been regarded as a complex phenomenon, which is generally characterized by depression of proteoglycan synthesis reduced water content, dysregulated extracellular matrix remodeling in the nucleus pulposus and inner annulus, as well as inappropriate collagen types (Sive et al., 2002; Le Maitre et al., 2005; Xu et al., 2016). Today, even though the knowledge of IDD has been greatly advanced, the outcomes of the treatment of IDD were unsatisfactory due to the complex mechanisms underlying IDD (Stokes and Iatridis, 2004; Wang et al., 2016; Zhang et al., 2017). Hence, it is still essential to elucidate the molecular mechanisms underlying the progression of IDD to develop novel and effective treatment of IDD.

Currently, genome-wide bioinformatics analysis as a popular approach for gaining insight into the overlapping and specific genetic mechanisms of disease is widely used. For example, Zhao et al. (2016) has identified 1854 long non-coding RNAs (lncRNAs) and 2804 protein-coding genes, which are differentially expressed in IDD group compared with the spinal cord injury group using RNA sequencing. Zhu et al. (2017) screened a total of 558 differentially expressed genes (DEGs), such as IL6, VEGFA, THBS1, and ITGA4, between the degenerated and normal samples based on array data of GSE42611. In our study, we used the gene expression profile of GSE70362 to further identify the molecular mechanism of IDD through identifying the DEGs, microRNAs (miRNAs), and interactions of miRNA-targets, as well as further analyzing their functions and pathways associated with the progression of IDD. The results obtained in this study may help to prove better insights into pathogenesis and define novel systemic biomarker for IDD.

2. Materials and Methods

2.1. Gene expression profile data

Microarray data sets with the access number GSE70362 were downloaded from the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo). The discs with Thompson degeneration grades I–II were assigned to the nondegenerative (control) group, while those having grades III–V were assigned to the degenerative group. The tissues were obtained from 24 donors, 8 of which (4 grade I, 2 grades I–II, and 2 grade II) were subject to the control group and the others (n = 16, including 6 grade III, 6 grade IV, and 4 grade V) were the degenerative group. The platform was GPL17810 (Affymetrix Human Genome U133 Plus 2.0 Array).

2.2. Data preprocessing and mining

The original CEL files were preprocessed by the robust multiarray average algorithm (Irizarry et al., 2003) from the “affy” package (Gautier et al., 2004). The preprocessing steps involve background adjustment, quantile normalization, finally summarization and probe ID to gene symbol, and after that the expression profile matrix was extracted. When different probes were mapping to the same gene symbol, the mean value of probes represents the expression value of gene. As a result, a total of 20,545 remaining genes were obtained.

2.3. Cluster analysis for DEGs

Noise robust soft clustering of gene expression data between control and degenerative groups was conducted by the Mfuzz package (Futschik and Carlisle, 2005; Kumar and Futschik, 2007). Then the fuzzy C-means clustering algorithm was used for gene ranking and clustering based on expression values. Multiple clustering results, including up- and downregulated DEGs, were obtained based on minSTD = 0.20, acore = 0.5, and cOptimal = 9.

2.4. Enrichment analysis and functional annotation

The Database for Annotation, Visualization, and Integrated Discovery clustering algorithms (Huang et al., 2007) were used for functional and pathway enrichment analysis of DEGs. The function and pathway enrichment analyses were conducted by the gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathways (Kanehisa and Goto, 2000). The threshold for the number of enriched genes was count ≥2.

2.5. Prediction of protein–protein interactions

Protein–protein interaction (PPI) between proteins was analyzed using the Search Tool for the Retrieval of Interacting Genes/Proteins (Mering et al., 2003) online tool (http://string-db.org), and the threshold for protein pairs was required confidence (combined score) >0.4.

2.6. Prediction of miRNA-target gene pair

Putative miRNAs that regulate genes in PPI network were predicted using Web-based Gene Set Analysis Tool Kit (http://bioinfo.vanderbilt.edu/webgestalt; Zhang et al., 2005). The threshold for the number of enriched genes was ≥4. The interactions of adjusted p < 0.05 based on Benjamini and Hochberg (Bader and Hogue, 2003) were selected, and the miRNA regulatory network was constructed.

2.7. Integration of PPI network and miRNA regulatory network

The PPI network and the miRNA regulatory network were integrated using the Cytoscape software (Shannon et al., 2003). Most biological networks have the common feature of scale-free property. The connectivity degree distribution of individual nodes (genes) in the network was analyzed to identify network hubs.

2.8. Network module detection

The MCODE Plugin (Bader and Hogue, 2003) of Cytoscape software (www.cytoscape.org) was implied to identify modules from the integrated network (score >5). In addition, GO and pathway enrichment analysis for genes in modules were performed.

3. Results

3.1. Gene clusters

A total of nine clusters were obtained, among which three clusters exhibiting an overall downward or upward trend which was coincident with the Thompson grades were further analyzed. As shown in Table 1, there were 127, 122, and 154 genes in clusters 2, 8, and 9, respectively. Moreover, the expression of genes in cluster 2 was significantly increased, while those in clusters 8 and 9 were significantly decreased, respectively (Fig. 1).

FIG. 1.

Three clusters obtained after genetic clustering. Every curve represents a gene. The darker color of the curve corresponds to a higher compliance with the overall gene expression trend.

Table 1.

Gene Clustering Analysis

Clusters	Tendency	Gene counts
Cluster 2	Up	127
Cluster 8	Down	122
Cluster 9	Down	154

3.2. Functional terms enriched by up- and downregulated genes

The results of functional enrichment are shown in Table 2. The upregulated genes were mainly associated with biological process (BP) of protein transport, establishment of protein localization, protein localization, and pathway of MAPK signaling pathway, while the downregulated genes were involved in GO-BP of membrane lipid metabolic process, regulation of endocytosis and positive regulation of endocytosis, and pathways of lysosome, drug metabolism, and fatty acid elongation in mitochondria.

Table 2.

The Functional Enrichment Results of Up- and Downregulated Differentially Expressed Genes

Category	Term	Description	Count	p
Upregulated genes
GOTERM_BP	GO:0015031	Protein transport	14	0.00170
GOTERM_BP	GO:0045184	Establishment of protein localization	14	0.00184
GOTERM_BP	GO:0008104	Protein localization	15	0.00220
GOTERM_BP	GO:0001666	Response to hypoxia	5	0.0127
GOTERM_BP	GO:0070482	Response to oxygen levels	5	0.0150
KEGG_PATHWAY	hsa04010	MAPK signaling pathway	6	0.0615
Downregulated genes
GOTERM_BP	GO:0006643	Membrane lipid metabolic process	5	0.0119
GOTERM_BP	GO:0030100	Regulation of endocytosis	4	0.0291
GOTERM_BP	GO:0045807	Positive regulation of endocytosis	3	0.0398
GOTERM_BP	GO:0009791	Postembryonic development	4	0.0458
GOTERM_BP	GO:0006123	Mitochondrial electron transport, cytochrome c to oxygen	2	0.0535
KEGG_PATHWAY	hsa04142	Lysosome	5	0.0370
KEGG_PATHWAY	hsa00983	Drug metabolism	3	0.0783
KEGG_PATHWAY	hsa00062	Fatty acid elongation in mitochondria	2	0.0834

BP, biological process; count, the number of genes enriched in the term; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

3.3. The prediction of miRNA–target interactions and integrated network construction

The significant miRNA–gene interactions are shown in Table 3, and among them, miR-124a, miR-145, miR-498, miR-506, and miR-96 had more targets. Next, the integrated miRNA–target gene regulatory network was constructed. As illustrated in Figure 2, the top 6 hub nodes of miRNAs in the network were miR-124a (degree = 18), miR-506 (degree = 18), miR-17-5p (degree = 14), miR-20a (degree = 14), miR-106a (degree = 14), and miR-20b (degree = 14). The top 6 hub genes were HSP90B1 (Heat Shock Protein 90 kDa Beta Family Member 1, degree = 18), FOS (FBJ Murine Osteosarcoma Viral Oncogene Homolog, degree = 17), IMPDH2 (Inosine 5′-monophosphate dehydrogenase 2, degree = 16), RPS4Y1 (Ribosomal Protein S4, Y-Linked 1, degree = 16), RPL15 (Ribosomal Protein L15, degree = 14), and UTY (Ubiquitously Transcribed Tetratricopeptide Repeat Containing, Y-Linked, degree = 14).

FIG. 2.

The integrated miRNA-target gene regulatory network. Red and green circles indicate up- and downregulated genes, respectively. The blue rhombus represents miRNA. T-shaped edges refer to the regulatory relationships between miRNAs and their targets; the others represent the interactions between proteins. The node size is in accordance with the degree. miRNA, microRNA.

Table 3.

The Significant MicroRNA–Target Gene Interactions

miRNA	Count	Adj. p value
hsa_TGCCTTA, miR-124a	18	0.0001
hsa_AACTGGA, miR-145	11	0.0003
hsa_GCTTGAA, miR-498	7	0.0016
hsa_GTGCCTT, miR-506	18	0.0016
hsa_GTGCCAA, miR-96	10	0.0065
hsa_GACTGTT, miR-212, miR-132	7	0.0074
hsa_GCACTTT, miR-17-5p, miR-20a, miR-106a, miR-106b, miR -20b, miR -519d	14	0.0081
hsa_ATGTCAC, miR-489	5	0.0103
hsa_CAGGGTC, miR-504	5	0.0103

p Value was adjusted by the multiple test adjustment (adj.p). Count indicates the number of genes regulated by miRNA.

miRNA, microRNA.

3.4. Modules obtained from the integrated network

With the threshold of score >5, a total of two modules were obtained from the integrated network (Fig. 3). There were 12 nodes and 37 interaction pairs in module 1 while the module 2 consisted of 7 nodes and 17 interaction pairs. Next, enrichment analysis of the genes in two modules is listed in Table 4. And the results showed that the cluster 1 was associated with signal peptide processing, cotranslational protein targeting to membrane, translation, protein targeting to membrane and peptide metabolic process, and pathway of Ribosome. Genes in cluster 2 were mainly enriched in GO-BP terms of protein folding.

FIG. 3.

Modules obtained from the integrated network. Red and green circles indicate up- and downregulated genes, respectively. Module 1 (A) and module 2 (B).

Table 4.

Gene Ontology Terms and Pathways Enriched by Genes in Two Modules

Category	Term	Description	Count	p
Cluster 1
GOTERM_BP_FAT	GO:0006465	Signal peptide processing	2	0.00531
GOTERM_BP_FAT	GO:0006613	Cotranslational protein targeting to membrane	2	0.00800
GOTERM_BP_FAT	GO:0006412	Translation	3	0.0192
GOTERM_BP_FAT	GO:0006612	Protein targeting to membrane	2	0.0217
GOTERM_BP_FAT	GO:0006518	Peptide metabolic process	2	0.0347
KEGG_PATHWAY	hsa03010	Ribosome	2	0.0505
Cluster 2
GOTERM_BP_FAT	GO:0006457	Protein folding	3	0.00177

4. Discussion

IDD is one of the common causes of low back pain, which affects most adults and is a major health problem in the world (Osborn and Smith, 1998). In this study, the gene expression profile generated from nondegenerative and degenerative intervertebral discs was analyzed. The results showed that the upregulated DEGs between two groups were associated with protein transport and localization, while the downregulated DEGs were enriched in membrane lipid metabolic process and endocytosis. After the integrated analysis of protein–protein and miRNA–target interactions, some hub genes such as HSP90B1, RPS4Y1, RPL15, and UTY, as well as hub miRNAs, including miR-124a and miR-506, were screened. Moreover, modules in the integrated network were functionally associated with protein targeting, peptide processing, and metabolic process, as well as protein folding. These genes and BPs may be involved in the progression of IDD, which could provide the theoretical basis to develop novel and effective treatment of IDD.

Currently, an increasing number of genes have been regarded to be involved in IDD development (Wu et al., 2015; Yan et al., 2015). In the present study, the DEGs, such as HSP90B1, RPS4Y1, RPL15, and UTY, between nondegenerative and degenerative intervertebral disc groups were identified as the key genes with the higher degrees. Actually, protein encoded by RPS4Y1 catalyzes protein synthesis, and increased protein synthesis can occur during disc degeneration (Cs-Szabo et al., 2002; Le Maitre et al., 2007), suggesting that RPS4Y1 plays an important role in IDD development. Similarly, the DEGs of HSP90B1 in cluster 2 were mainly found to be associated with protein folding, of which the encoded protein functions in the stabilizing and folding of other proteins and involves in a variety of pathogenic states (Bloor et al., 2010). In addition, an increasing number of abnormally expressed miRNAs, such as miR-193a-3p and miR-27b, serve important roles in the regulation of IDD through targeting matrix metalloproteinase (MMP)14 and MMP13, respectively (Ji et al., 2016; Li et al., 2016). In accordance with previous studies, some hub miRNAs, including miR-124a and miR-506, in our study could also be involved in the mechanism of IDD. However, little correlation between miRNAs and hub genes was identified. Thus the miRNA expression profile may be needed to integrate differential miRNA expression in gene expression in our next study.

Previous studies suggested that occlusion of vertebral end plate openings limits the transport of nutrients and metabolites and affects the loss and synthesis of matrix structural proteins, resulting in disc degeneration (Bibby et al., 2002; Benneker et al., 2005). Moreover, the site-specific localization of cathepsins was involved in the pathomechanism of IDD (Ariga et al., 2001), and the MAPK pathway plays a role in osmoregulation in intervertebral disc cells through TonEBP and its target genes (Tsai et al., 2007). The receptor-mediated endocytosis is mediated by the MAPK signaling pathway (Gong et al., 2007). In our study, the gene cluster analysis revealed that many genes, especially those in clusters 2, 8, and 9, were gradually changed in gene expression and were coincident with the degenerative grade. In the following analysis, upregulated DEGs that may be involved in IDD were enriched in protein transport and localization and MAPK signaling pathway. In addition, membrane lipid metabolic process and endocytosis enriched by downregulated DEGs may also be important for the development of IDD. Actually, even though several studies have proposed that serum lipid levels may be a risk factor for intervertebral disc pathology (Longo et al., 2011), and lipid peroxidation contributes to degeneration with aging (Ames et al., 1993), there is no direct evidence about the correlation between membrane lipid metabolism and IDD. Hence, our data could provide the theoretical basis for the related study in future.

In conclusion, protein translation, folding, transport, targeting, and localization could be greatly dysregulated during IDD. Some genes, including RPS4Y1, UTY, RPL15, and HSP90B1, were abnormally expressed and may be important for the development of this disease. However, there are some limitations in the current analysis that should be of concern, including the data generated from a single platform. Moreover, some validation experiments are still needed to confirm the results.

Footnotes

Author Disclosure Statement

The authors declare they have no competing financial interests.

Funding Information

No funding was received for this article.

References

Ames

B.N.

, Shigenaga

M.K.

, and Hagen

T.M.

1993. Oxidants, antioxidants, and the degenerative diseases of aging. Proc. Natl. Acad. Sci. USA. 90, 7915–7922.

Ariga

, Yonenobu

, Nakase

, et al. 2001. Localization of cathepsins D, K, and L in degenerated human intervertebral discs. Spine. 26, 2666–2672.

Bader

G.D.

, and Hogue

C.W.

2003. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinformatics. 4, 2.

Benneker

L.M.

, Heini

P.F.

, Alini

, et al. 2005. 2004 Young Investigator Award Winner: Vertebral endplate marrow contact channel occlusions and intervertebral disc degeneration. Spine. 30, 167–173.

Bibby

S.R.

, Fairbank

J.C.

, Urban

M.R.

, et al. 2002. Cell viability in scoliotic discs in relation to disc deformity and nutrient levels. Spine. 27, 2220–2227.

Bloor

, Maelfait

, Krumbach

, et al. 2010. Endoplasmic reticulum chaperone gp96 is essential for infection with vesicular stomatitis virus. Proc. Natl. Acad. Sci. USA. 107, 6970–6975.

Cs-Szabo

, Ragasa-San Juan

, Turumella

, et al. 2002. Changes in mRNA and protein levels of proteoglycans of the anulus fibrosus and nucleus pulposus during intervertebral disc degeneration. Spine. 27, 2212–2219.

Futschik

M.E.

, and Carlisle

2005. Noise-robust soft clustering of gene expression time-course data. J. Bioinform. Comput. Biol. 3, 965–988.

Gautier

, Cope

, Bolstad

B.M.

, et al. 2004. affy—analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 20, 307–315.

10.

Gong

, Wei

, Li

, et al. 2007. p38 MAPK signaling pathway is involved in the regulation of receptor mediated endocytosis. 23, P1388–P1391.

11.

Huang

D.W.

, Sherman

B.T.

, Tan

, et al. 2007. The DAVID Gene Functional Classification Tool: A novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 8, R183.

12.

Irizarry

R.A.

, Hobbs

, Collin

, et al. 2003. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 4, 249–264.

13.

Jensen

C.E.

, Riis

, Petersen

K.D.

, et al. 2017. Economic evaluation of an implementation strategy for the management of low back pain in general practice. Pain. 158, 891–899.

14.

M.-L.

, Zhang

X.-J.

, Shi

P.-L.

, et al. 2016. Downregulation of microRNA-193a-3p is involved in invertebral disc degeneration by targeting MMP14. J. Mol. Med. 94, 457–468.

15.

Kanehisa

, and Goto

2000. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res. 28, 27–30.

16.

Kumar

, and Futschik

M.E.

2007. Mfuzz: A software package for soft clustering of microarray data. Bioinformation. 2, 5.

17.

Le Maitre

, Pockert

, Buttle

, et al. 2007. Matrix synthesis and degradation in human intervertebral disc degeneration. Biochem. Soc. Trans. 35, 652–655.

18.

Le Maitre

C.L.

, Freemont

A.J.

, and Hoyland

J.A.

2005. The role of interleukin-1 in the pathogenesis of human intervertebral disc degeneration. Arthritis Res. Ther. 7, R732.

19.

H.-R.

, Cui

, Dong

Z.-Y.

, et al. 2016. Downregulation of miR-27b is involved in loss of type II collagen by directly targeting matrix metalloproteinase 13 (MMP13) in human intervertebral disc degeneration. Spine. 41, E116–E123.

20.

Longo

U.G.

, Denaro

, Spiezia

, et al. 2011. Symptomatic disc herniation and serum lipid levels. Eur. Spine J. 20, 1658–1662.

21.

Mering

C.V.

, Huynen

, Jaeggi

, et al. 2003. STRING: A database of predicted functional associations between proteins. Nucleic Acids Res. 31, 258–261.

22.

Millecamps

, Tajerian

, Naso

, et al. 2012. Lumbar intervertebral disc degeneration associated with axial and radiating low back pain in ageing SPARC-null mice. Pain. 153, 1167–1179.

23.

Osborn

, and Smith

J.A.

1998. The personal experience of chronic benign lower back pain: An interpretative phenomenological analysis. Br. J. Health Psychol. 3, 65–83.

24.

Raj

P.P.

2008. Intervertebral disc: Anatomy-physiology-pathophysiology-treatment. Pain Pract. 8, 18–44.

25.

Shannon

, Markiel

, Ozier

, et al. 2003. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504.

26.

Sive

, Baird

, Jeziorsk

, et al. 2002. Expression of chondrocyte markers by cells of normal and degenerate intervertebral discs. Mol. Pathol. 55, 91.

27.

Stokes

I.A.

, and Iatridis

J.C.

2004. Mechanical conditions that accelerate intervertebral disc degeneration: Overload versus immobilization. Spine. 29, 2724–2732.

28.

Tsai

T.T.

, Guttapalli

, Agrawal

, et al. 2007. MEK/ERK signaling controls osmoregulation of nucleus pulposus cells of the intervertebral disc by transactivation of TonEBP/OREBP. J. Bone Miner. Res. 22, 965–974.

29.

Wang

, Wang

, Yan

, et al. 2016. MiR-138-5p promotes TNF-α-induced apoptosis in human intervertebral disc degeneration by targeting SIRT1 through PTEN/PI3K/Akt signaling. Exp. Cell Res. 345, 199–205.

30.

, Wang

, Chen

, et al. 2015. Resveratrol could reverse the expression of SIRT1 and MMP-1 in vitro. Genet. Mol. Res. 14, 12386–12393.

31.

J., E, X.Q.

, Wang

N.X.

, et al. 2016. BMP 7 enhances the effect of BMSC s on extracellular matrix remodeling in a rabbit model of intervertebral disc degeneration. FEBS J. 283, 1689–1700.

32.

Yan

, Yu

, Zhang

, et al. 2015. Lumbar disc degeneration is facilitated by MiR-100-mediated FGFR3 suppression. Cell. Physiol. Biochem. 36, 2229–2236.

33.

Zhang

, Kirov

, and Snoddy

2005. WebGestalt: An integrated system for exploring gene sets in various biological contexts. Nucleic Acids Res. 33, W741–W748.

34.

Zhang

, Zhao

, Shen

, et al. 2016. Molecular mechanisms of cell death in intervertebral disc degeneration. Int. J. Mol. Med. 37, 1439–1448.

35.

Zhang

W.L.

, Chen

Y.F.

, Meng

H.Z.

, et al. 2017. Role of miR-155 in the regulation of MMP-16 expression in intervertebral disc degeneration. J. Orthop. Res. 35, 1323–1334.

36.

Zhao

, Lu

, Wang

, et al. 2016. Genome-wide identification of long noncoding RNAs in human intervertebral disc degeneration by RNA sequencing. Biomed Res. Int. 2016, 3684875.

37.

Zhu

, Chen

, Jiao

, et al. 2017. Identification of critical genes in nucleus pulposus cells isolated from degenerated intervertebral discs using bioinformatics analysis. Mol. Med. Rep. 16, 553–564.