Gene expression profile analysis of human intervertebral disc degeneration

GENOMICS AND BIOINFORMATICS

RESEARCH ARTICLE

Gene expression profile analysis of human intervertebral disc degeneration

Kai Chen^*; Dajiang Wu^*; Xiaodong Zhu; Haijian Ni; Xianzhao Wei; Ningfang Mao; Yang Xie; Yunfei Niu; Ming Li

Department of Orthopedics, Changhai Hospital affiliated to the Second Military Medical University, Shanghai, China

^{Send correspondence to}

ABSTRACT

In this study, we used microarray analysis to investigate the biogenesis and progression of intervertebral disc degeneration. The gene expression profiles of 37 disc tissue samples obtained from patients with herniated discs and degenerative disc disease collected by the National Cancer Institute Cooperative Tissue Network were analyzed. Differentially expressed genes between more and less degenerated discs were identified by significant analysis of microarray. A total of 555 genes were significantly overexpressed in more degenerated discs with a false discovery rate of < 3%. Functional annotation showed that these genes were significantly associated with membrane-bound vesicles, calcium ion binding and extracellular matrix. Protein-protein interaction analysis showed that these genes, including previously reported genes such as fibronectin, COL2A1 and f-catenin, may play key roles in disc degeneration. Unsupervised clustering indicated that the widely used morphology-based Thompson grading system was only marginally associated with the molecular classification of intervertebral disc degeneration. These findings indicate that detailed, systematic gene analysis may be a useful way of studying the biology of intervertebral disc degeneration.

Keywords: genes, intervertebral disc degeneration, molecular classification, protein-protein interaction.

Introduction

The intervertebral disc (IVD) is a cartilaginous structure that shows degenerative and age-related changes earlier than other connective tissues in the body (Urban and Roberts, 2003). Disc degeneration occurs when its cells die or become dysfunctional, especially in an acidic environment. During degeneration, the IVD becomes dehydrated and vascularized, and there is an ingrowth of nerves (Hughes et al., 2012). As a multi-factorial process, disc degeneration is influenced by factors such as genetic predisposition, lifestyle conditions, including obesity, occupation, smoking and alcohol consumption (Samartzis et al., 2012), aging phenomena (Gruber and Hanley, 2007) and other health factors, such as diabetes. Disc degeneration is important clinically because of a strong association between disc degeneration and back pain (Luoma et al., 2000), which is a major public health problem in modern industrialized societies. In addition, the altered physiology of the IVD during degeneration is associated with many other clinical symptoms or conditions, such as lower limb pain, paraesthesia, spinal stenosis and disc herniation (Hughes et al., 2012).

Previous studies have identified some molecular markers associated with disc degeneration. For example, asporin was found to be expressed at higher levels in degenerated human IVD (Gruber et al., 2009). The mRNA expression levels of ADAMTS-7 and ADAMTS-12 were significantly higher in endplate cells of degenerative discs compared with those in non-degenerative discs (Zhang et al., 2012). In addition, genetic polymorphisms in collagen I, collagen IX, collagen XI, aggrecan, extracellular matrix-degrading enzymes, inflammatory cytokines, such as

IL (interleukin)-1, IL-6 and TNFa (tumor necrosis factors a), Fas/FasL and vitamin D receptors have been associated with the development of intervertebral disc degeneration (IVDD) (Kalb et al., 2012).

In this study, we compared the gene expression profile data of IVDD patients with that of healthy controls. The gene profiles were downloaded from the gene expression omnibus (GEO) database and the intervertebral disc data were classified into two groups according to the Thompson grading system (Thompson et al., 1990).

Materials and Methods

Microarray data analysis

Microarray datasets were downloaded from the GEO database (accession numbers GSE15227 (Gruber et al., 2009) and GSE23130 (Gruber et al., 2012)). The signal intensities were re-calculated with the robust multi-array average (RMA) (Irizarry et al., 2003) using the Entrez Gene-based center custom Chip Description File (CDF) developed by University of Michigan (Dai et al., 2005). The single sample classified as Thompson grade I was considered as healthy (Gruber et al., 2009) and ignored before merging the GSE15227 and GSE23130 data. ComBat (Johnson et al., 2007) was applied to the merged dataset to remove potential unwanted batch effects. The samples were divided into two groups based on the Thompson grading system. Grade IV and V samples were considered as more degenerated (Gruber et al., 2009). Significant analysis of microarray (SAM) (Tusher et al., 2001) was applied to identify differentially expressed genes using a delta of 1. Functional annotation was done on DAVID (Huang da et al., 2009) using significantly differentially expressed genes.

Protein-protein interaction analysis

Differentially expressed genes were submitted to STRING 9.0 (Szklarczyk et al., 2011) to detect possible protein-protein interactions. EntrezGene IDs were translated into official gene symbols by DAVID before submitting. A threshold score of 0.4 was applied as default. Only those interactions that were proven experimentally were considered. A simple Perl script was used to exclude genes/proteins not in the query gene list.

Cluster analysis

The data for each gene was median-centered before clustering. Centered correlation was applied as a similarity metric and average linkage as a clustering method. The analysis was done in Cluster 3.0 and a heatmap was generated by TreeView. A one-sided Fisher's exact test was used to assess the consistence of the molecular and Thompson grade classifications.

Results

Differentially expressed genes in IVD degeneration

Eleven disc samples (eight Grade IV and three Grade V) were classified as the more degenerated group while the remaining 26 samples (10 Grade II and 16 Grade III) were classified as the less degenerated group. One thousand sample permutations were used to estimate the false discovery rate (FDR). When delta was set to 1,555 genes were overexpressed in the more degenerated disc samples (Table 1) and 98 genes were underexpressed in the more degenerated disc samples (Table 2). Both groups had an FDR < 3%. The overexpressed gene list was referred to as IVDD_UP and the underexpressed gene list as IVDD_DN. Several of the IVDD_UP genes have been reported before. For example, asporin, which is present in the cartilage extracellular matrix (ECM) (Henry et al., 2001) and may be genetically associated with osteoarthritis (Loughlin, 2005), showed higher expression levels in Thompson Grade IV human discs compared to lower grade discs (Gruber et al., 2009). A highly conserved secreted serine protease, HTRA1, which degrades numerous extracellular matrix proteins, was also found in IVDD_UP. HTRA1 mRNA and protein levels are significantly elevated in degenerated disc tissue (Tiaden et al., 2012). The gene expression of cathepsin K, a cysteine protease that cleaves the triple helical domains of collagen types I and II, was significantly greater in more degenerated discs (grades III and IV) compared to healthier discs (grades I and II) (Gruber et al., 2011).

Functional annotation of differentially expressed genes

Differentially expressed genes were submitted to DAVID (Huang da et al., 2009) for functional annotation. Databases for disease information, functional categories, gene ontology, curated pathways and protein domains were included to provide a comprehensive knowledge source. For the IVDD_UP genes, 195 functional annotation clusters were identified. The top three clusters ranked by the enrichment score were mainly associated with membranebounded vesicles, calcium ion binding and extracellular matrix, with enrichment scores of 4.7, 3.61 and 3.25, respectively. Noticeably, these enriched genes were also related to the biological process of skeletal system development, including cartilage development. Among the IVDD_DN genes, there were only 19 clusters with the highest enrichment score being only 0.87; this score was for a cluster involving mainly zinc fingers and the Krueppelassociated box.

Protein-protein interactions among differentially expressed genes

STRING showed that 159 experimentally proven protein-protein interactions were formed by 175 genes in IVDD_UP (Figure 1). Fibronectin 1, COL2A1 (collagen, type II, α1) and β-catenin were the three genes/proteins that had the most interacting partners (Figure 2). Of these 175 genes, 107 had only one interacting partner in IVDD_UP and thus resembled the phenomenon of scale-free biological networks (Barabasi and Albert, 1999). Surprisingly, no interactions were observed amongst IVDD_DN genes.

Cluster analysis of IVD tissues

Unsupervised hierarchical clustering was done on 22 disc tissues of GSE23130. The experimental results obtained using all 18,818 genes or only 1,984 genes with a minimum sample standard deviation of 1 were highly similar (Figure 3). A one-tailed Fisher's exact test showed a marginally significant association between the Thompsonbased classification and the unsupervised molecular classification (p = 0.04799 for the 18,818-gene case and 0.09133 for the 1,984-gene case).

Discussion

Our microarray analysis of the gene expression profile of degenerated IVD detected hundreds of differentially expressed genes that may be associated with IVD degeneration. Functional annotation supported the biological relevance of our findings. Protein-protein interaction analysis of IVDD_UP revealed genes that may have key roles in IVD degeneration. Fibronectin 1 (FN1), a molecular marker for fibrosis (Leask and Abraham, 2004), was upregulated in punctured mouse tail discs that showed progressive degenerative changes and were fibrocartilaginous (Yang et al., 2009). Earlier studies reported that the degradation fragments of fibronectin (Fn-f) accumulate in the disc during degeneration (Oegema et al., 2000) and induce IVD degeneration in rabbits (Greg Anderson et al., 2003). Given that HTRA1 can induce fibronectin proteolysis and one of the resultant fragments was found to be increased in HTRA1-treated IVD cell cultures as well as in disc tissue from patients with IVD degeneration, it was hypothesized that HTRA1 promoted IVD degeneration through the proteolytic cleavage of fibronectin and subsequent activation of resident disc cells (Tiaden et al., 2012). The overexpression of these two genes was consistent with this hypothesis. Interestingly, fibronectin mRNA and protein levels were significantly down-regulated following the upregulation of HTRA1 in rhesus monkey choroid-retina endothelial cells (RF/6A) and human umbilical vein endothelial cells (HUVECs) (Jiang et al., 2012), suggesting that some other unknown factors may be involved in the interaction between HTRA1 and fibronectin in IVD tissues.

The IVDD_UP group of genes contained four collagen genes: COL11A1, COL1A2, COL3A1 and COL2A1. COL2A1 interacted with eight other proteins encoded by genes in IVDD_UP. Mutations in COL2A1 give rise to a spectrum of phenotypes that predominantly affecting cartilage and bone. These chondrodysplasias are typically characterized by a disproportionately short stature, eye abnormalities, cleft palate and hearing loss. Increasing evidence has also implicated COL2A1 mutations in individuals with isolated degenerative joint disease (Kannu et al., 2010). The relationship between COL2A1 gene polymorphisms and knee osteoarthritis were also investigated in Han Chinese women (Xu et al., 2011).

Mutations in the COL11A1 gene (also known as STL2) have been identified in some people with Stickler syndrome (Martin et al., 1999). Other mutations in this gene cause segments of DNA to be skipped when the protein is being made, resulting in an abnormally short proalpha chain. Alterations in COL11A1 impair collagen function and can lead to hearing loss, tearing of the lining of the eye, and bone and joint abnormalities (Rodriguez et al., 2004). COL1A2 is a protein found in most connective tissues. Mutations in this gene are associated with asteogenesis imperfecta, Ehlers-Danlos syndrome, idiopathic osteoporosis and atypical Marfan syndrome (Vasan et al., 1991; Ward et al., 2001; Hartikka et al., 2004). However, the symptoms associated with mutations in this gene tend to be less severe than with mutations in the gene for a1 type I collagen since the α2 form is less abundant (Bou-Gharios et al., 2004). Multiple messages for this gene result from multiple polyadenylation signals, a feature shared with most of the other collagen genes (Zuo et al., 2012).

COL3A1 is a protein that, in humans, is encoded by the COL3A1 gene located on chromosome 2 (Cutting et al., 1990). COL3A1 is a fibrillar collagen found in extensible connective tissues such as skin, lung and the vascular system, frequently in association with type 1 collagen. Mutations in this gene often lead to the exclusion of multi-exons (Superti-Furga et al., 1988). Other studies have shown that (1) IVD cells respond strongly to changes in the osmotic environment by altering their mRNA expression, e.g.,human cells cultured over five days showed increased expression of aggrecan and collagen II in the nucleus and annulus cells in the presence of elevated osmolarity (Wuertz et al., 2007), and (2) there is an increase in collagen-II, aggrecan and Sox-9 protein expression in the nucleus pulposus (NP) and anulus fibrosus (AF) regions of discs from running exercised rats compared with non-exercised controls (Brisby et al., 2010). Our result may indicate an as yet unknown role for COL2A1 in degenerated IVD.

Proper regulation of Wnt/β-catenin signaling is required for development and organization of the IVD (Kondo et al., 2011). Furthermore, β-catenin was overexpressed in IVD extracted from patients with IVD degeneration and X-ray and micro-CT analyses revealed osteophyte formation and narrowing of the disc space in three-month-old β-catenin conditional activation (cAct) mice (Tang et al., 2012). Our results may provide clues for studying the biological mechanism of f-catenin in human disc tissue.

Unsupervised clustering has provided insights into the molecular heterogeneity of complex diseases such as cancer and is useful in cancer diagnosis and therapy. The application of this analysis in our study showed that the widely used morphology-based Thompson grading system was only marginally associated with the molecular classification of IVDD. This interesting finding indicates that there is scope for improving the clinical assessment of the progress of IVDD.

Received: January 11, 2013; Accepted: June 11, 2013.

• Barabasi AL and Albert R (1999) Emergence of scaling in random networks. Science 286:509-512.
• Bou-Gharios G, Ponticos M, Rajkumar V and Abraham D (2004) Extra-cellular matrix in vascular networks. Cell Prolif 37:207-220.
• Brisby H, Wei AQ, Molloy T, Chung SA, Murrell GA and Diwan AD (2010) The effect of running exercise on intervertebral disc extracellular matrix production in a rat model. Spine (Phila Pa 1976) 35:1429-1436.
• Cutting GR, McGinniss MJ, Kasch LM, Tsipouras P and Antonarakis SE (1990) Physical mapping by PFGE localizes the COL3A1 and COL5A2 genes to a 35-kb region on human chromosome 2. Genomics 8:407-410.
• Dai M, Wang P, Boyd AD, Kostov G, Athey B, Jones EG, Bunney WE, Myers RM, Speed TP, Akil H, et al. (2005) Evolving gene/transcript definitions significantly alter the interpretation of GeneChip data. Nucleic Acids Res 33:e175.
• Greg Anderson D, Li X, Tannoury T, Beck G and Balian G (2003) A fibronectin fragment stimulates intervertebral disc degeneration in vivo. Spine (Phila Pa 1976) 28:2338-2345.
• Gruber HE and Hanley EN (2007) Do we need biomarkers for disc degeneration? Biomark Insights 1:131-133.
• Gruber HE, Hoelscher G, Loeffler B, Chow Y, Ingram JA, Halligan W and Hanley Jr EN (2009a) Prostaglandin E1 and misoprostol increase epidermal growth factor production in 3D-cultured human annulus cells. Spine J 9:760-766.
• Gruber HE, Ingram JA, Hoelscher GL, Zinchenko N, Hanley r. EN and Sun Y (2009b) Asporin, a susceptibility gene in osteoarthritis, is expressed at higher levels in the more degenerate human intervertebral disc. Arthritis Res Ther 11:R47.
• Gruber HE, Ingram JA, Hoelscher GL, Zinchenko N, Norton HJ and Hanley Jr EN (2011) Constitutive expression of cathepsin K in the human intervertebral disc: New insight into disc extracellular matrix remodeling via cathepsin K and receptor activator of nuclear factor-kappaB ligand. Arthritis Res Ther 13:R140.
• Gruber HE, Hoelscher GL, Ingram JA and Hanley Jr. EN (2012) Genome-wide analysis of pain-, nerve-and neurotrophinrelated gene expression in the degenerating human annulus. Mol Pain 8:63.
• Hartikka H, Kuurila K, Körkkö J, Kaitila I, Grénman R, Pynnönen S, Hyland JC and Ala-Kokko L (2004) Lack of correlation between the type of COL1A1 or COL1A2 mutation and hearing loss in osteogenesis imperfecta patients. Hum Mutat 24:147-154.
• Henry SP, Takanosu M, Boyd TC, Mayne PM, Eberspaecher H, Zhou W, de Crombrugghe B, Hook M and Mayne R (2001) Expression pattern and gene characterization of asporin, a newly discovered member of the leucine-rich repeat protein family. J Biol Chem 276:12212-12221.
• Huang da W, Sherman BT and Lempicki RA (2009) Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 4:44-57.
• Hughes SP, Freemont AJ, Hukins DW, McGregor AH and Roberts S (2012) The pathogenesis of degeneration of the intervertebral disc and emerging therapies in the management of back pain. J Bone Joint Surg Br 94:1298-1304.
• Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U and Speed TP (2003) Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 4:249-264.
• Jiang J, Huang L, Yu W, Wu X, Zhou P and Li X (2012) Overexpression of HTRA1 leads to down-regulation of fibronectin and functional changes in RF/6A cells and HUVECs. PLoS One 7:e46115.
• Johnson WE, Li C and Rabinovic A (2007) Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8:118-127.
• Kalb S, Martirosyan NL, Kalani MY, Broc GG and Theodore N (2012) Genetics of the degenerated intervertebral disc. World Neurosurg 77:491-501.
• Kannu P, Bateman JF, Randle S, Cowie S, du Sart D, McGrath S, Edwards M and Savarirayan R (2010) Premature arthritis is a distinct type II collagen phenotype. Arthritis Rheum 62:1421-1430.
• Kondo N, Yuasa T, Shimono K, Tung W, Okabe T, Yasuhara R, Pacifici M, Zhang Y, Iwamoto M and Enomoto-Iwamoto M (2011) Intervertebral disc development is regulated by Wnt/beta-catenin signaling. Spine (Phila Pa 1976) 36:E513-E518.
• Leask A and Abraham DJ (2004) TGF-beta signaling and the fibrotic response. FASEB J 18:816-827.
• Loughlin J (2005) Polymorphism in signal transduction is a major route through which osteoarthritis susceptibility is acting. Curr Opin Rheumatol 17:629-633.
• Luoma K, Riihimaki H, Luukkonen R, Raininko R, Viikari-Juntura E and Lamminen A (2000) Low back pain in relation to lumbar disc degeneration. Spine (Phila Pa 1976) 25:487-492.
• Martin S, Richards AJ, Yates J, Scott JD, Pope M and Snead MP (1999) Stickler syndrome: further mutations in COL11A1 and evidence for additional locus heterogeneity. Eur J Hum Genet 7:807-814
• Oegema Jr. TR, Johnson SL, Aguiar DJ and Ogilvie JW (2000) Fibronectin and its fragments increase with degeneration in the human intervertebral disc. Spine (Phila Pa 1976) 25:2742-2747.
• Rodriguez R, Seegmiller R, Stark M and Bridgewater L (2004) A type XI collagen mutation leads to increased degradation of type II collagen in articular cartilage. Osteoarthritis Cartilage 12:314-320.
• Samartzis D, Karppinen J, Chan D, Luk KD and Cheung KM (2012) The association of lumbar intervertebral disc degeneration on magnetic resonance imaging with body mass index in overweight and obese adults: a population-based study. Arthritis Rheum 64:1488-1496.
• Superti-Furga A, Gugler E, Gitzelmann R and Steinmann B (1988) Ehlers-Danlos syndrome type IV: a multi-exon deletion in one of the two COL3A1 alleles affecting structure, stability, and processing of type III procollagen. J Biol Chem 263:6226-6232.
• Szklarczyk D, Franceschini A, Kuhn M, Simonovic M, Roth A, Minguez P, Doerks T, Stark M, Muller J, Bork P, et al. (2011) The STRING database in 2011: Functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res 39:D561-D568.
• Tang D, Xu G, Hou W, Ye X, Shi Q, Chen D and Wang Y (2012) β-Catenin plays a key role in lumbar intervertebral disc degeneration. Spine (Phila Pa 1976)18-37
• Thompson JP, Pearce RH, Schechter MT, Adams ME, Tsang IK and Bishop PB (1990) Preliminary evaluation of a scheme for grading the gross morphology of the human intervertebral disc. Spine (Phila Pa 1976) 15:411-415.
• Tiaden AN, Klawitter M, Lux V, Mirsaidi A, Bahrenberg G, Glanz S, Quero L, Liebscher T, Wuertz K, Ehrmann M, et al. (2012) Detrimental role for human high temperature requirement serine protease A1 (HTRA1) in the pathogenesis of intervertebral disc (IVD) degeneration. J Biol Chem 287:21335-21345.
• Tusher VG, Tibshirani R and Chu G (2001) Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci USA 98:5116-5121.
• Urban JP and Roberts S (2003) Degeneration of the intervertebral disc. Arthritis Res Ther 5:120-130.
• Vasan NS, Kuivaniemi H, Vogel BE, Minor R, Wootton J, Tromp G, Weksberg R and Prockop D (1991) A mutation in the pro alpha 2 (I) gene (COL1A2) for type I procollagen in Ehlers-Danlos syndrome type VII: Evidence suggesting that skipping of exon 6 in RNA splicing may be a common cause of the phenotype. Am J Hum Genet 48:305-317.
• Ward L, Lalic L, Roughley P and Glorieux F (2001) Thirty-three novel COL1A1 and COL1A2 mutations in patients with osteogenesis imperfecta types I-IV. Hum Mutat 17:434-434.
• Wuertz K, Urban JP, Klasen J, Ignatius A, Wilke HJ, Claes L and Neidlinger-Wilke C (2007) Influence of extracellular osmolarity and mechanical stimulation on gene expression of intervertebral disc cells. J Orthop Res 25:1513-1522.
• Xu P, Yao J and Hou W (2011) Relationships between COL2A1 gene polymorphisms and knee osteoarthritis in Han Chinese women. Mol Biol Rep 38:2377-2381.
• Yang F, Leung VY, Luk KD, Chan D and Cheung KM (2009) Injury-induced sequential transformation of notochordal nucleus pulposus to chondrogenic and fibrocartilaginous phenotype in the mouse. J Pathol 218:113-121.
• Zhang Q, Huang M, Wang X, Xu X, Ni M and Wang Y (2012) Negative effects of ADAMTS-7 and ADAMTS-12 on endplate cartilage differentiation. J Orthop Res 30:12381243.
• Zuo C, Wen F, Li M, Zhang X, Chen H, Wu K and Zeng R (2012) COL1A2 polymorphic markers confer an increased risk of neovascular age-related macular degeneration in a Han Chinese population. Mol Vision 18:1787-1793.

*

These authors contributed equally to this work.

Associate Editor: Houtan Noushmer

License information: This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Publication Dates

• Publication in this collection
  23 Aug 2013
• Date of issue
  2013

History

• Received
  11 Jan 2013
• Accepted
  11 June 2013