Gene expression profile analysis of human intervertebral disc degeneration.

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 β-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.

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 ). 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 ), 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 ), 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 ).

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 ). 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 ). 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 TNFα (tumor necrosis factors α), Fas/FasL and vitamin D receptors have been associated with the development of intervertebral disc degeneration (IVDD) (Kalb ).

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 ).

Materials and Methods

Microarray data analysis

Microarray datasets were downloaded from the GEO database (accession numbers GSE15227 (Gruber ) and GSE23130 (Gruber )). The signal intensities were re-calculated with the robust multi-array average (RMA) (Irizarry ) using the Entrez Gene-based center custom Chip Description File (CDF) developed by University of Michigan (Dai ). The single sample classified as Thompson grade I was considered as healthy (Gruber ) and ignored before merging the GSE15227 and GSE23130 data. ComBat (Johnson ) 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 ). Significant analysis of microarray (SAM) (Tusher ) was applied to identify differentially expressed genes using a delta of 1. Functional annotation was done on DAVID (Huang da ) using significantly differentially expressed genes.

Protein-protein interaction analysis

Differentially expressed genes were submitted to STRING 9.0 (Szklarczyk ) 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 over-expressed 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 ) 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 ). 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 ). 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 ).

Table 1

Top 30 IVDD_UP genes with p < 0.05 and FDR < 3%.

Gene name	Score	Fold change	q-value
BCAS3	5.525645110	2.015911968	0
AP2S1	4.993647193	2.695484514	0
RTN3	4.831739884	1.906392737	0
MAT2B	4.795801549	2.661280400	0
SMOC2	4.722070979	3.529131297	0
TMEM214	4.609155146	2.397158886	0
SETD3	4.553095283	1.580288009	0
RCAN1	4.515015896	2.901977015	0
NUDT15	4.387565829	1.611085005	0
PDGFRA	4.352020487	2.013095239	0
HOXC6	4.312467370	4.046628242	0
DAXX	4.270628387	1.822846390	0
SLC40A1	4.255751181	5.142171070	0
RYFY1	4.250428720	2.629236672	0
IFITM1	4.237500892	2.139999154	0
CHMP1A	4.220758109	2.250578926	0
SF4	4.213669108	1.548679838	0
GJA1	4.184343713	5.752785110	0
ZFYVE26	4.095399931	1.536943270	0
CHT3L1	4.080014011	5.398762712	0
C3ORF1	4.069293869	2.752723264	0
CPN2	4.067581072	4.689686574	0
POLR2G	4.050051549	2.546570992	0
LGALS1	4.027157407	4.896502415	0
ITM2B	4.026549335	3.916969611	0
PTTG1	4.006722096	4.387834233	0
CDK2AP1	3.968504177	6.601587442	0
TMEM9	3.967576485	1.540136100	0
CLU	3.947513637	7.942443254	0
RHOC	3.942703220	3.185787316	0

Table 2

Top 30 IVDD_DN genes with p < 0.05 and FDR < 3%.

Gene name	Score	Fold change	q-value
KRTAP4-11	−4.465760325	0.478119384	0
OR5E1P	−4.306673904	0.598139097	0
TTTY14	−4.266771085	0.697225016	0
PPY	−4.230640542	0.327213756	0
C100RF25	−4.209818145	0.58437279	0
DPCR1	−4.101123182	0.700400802	0
LOC1002893	−4.052210058	0.560921078	0
DNAH8	−3.978399529	0.708059442	0
CSMD2	−3.948145754	0.685029482	0
KRTAP4-6	−3.870547452	0.669765942	0
HOXD11	−3.810598074	0.706993483	0.417452925
TMEM102	−3.798006024	0.71228122	0.417452925
KRTAP4-8	−3.76590512	0.625274173	0.417452925
KRTAP4-4	−3.745342914	0.733361857	0.417452925
ZNF440	−3.7332524	0.644960742	0.417452925
OR2M4	−3.721057881	0.656210145	0.417452925
TKTL1	−3.659388232	0.704826703	0.538764031
KRTAP-7	−3.657243192	0.692592119	0.538764031
KRT76	−3.633206767	0.533949119	0.538764031
C80RF12	−3.629644489	0.677570179	0.538764031
LOC1001312	−3.562150167	0.721969038	0.538764031
TFE3	−3.55931541	0.789925298	0.538764031
FRS3	−3.540271969	0.622015836	0.538764031
HHIP	−3.525281501	0.717589029	0.538764031
HTR3C	−3.522633996	0.603942904	0.538764031
ELAC1	−3.519001164	0.785612099	0.538764031
LOC401097	−3.492245688	0.766487773	1.000561772
TRAV24	−3.487013115	0.656667894	1.000561772
FMO1	−3.457758868	0.658106046	1.000561772
WDR53	−3.456389239	0.626137741	1.000561772

Functional annotation of differentially expressed genes

Differentially expressed genes were submitted to DAVID (Huang da ) 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 membrane-bounded 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 Krueppel-associated 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.

Figure 1

Protein-protein interaction network for IVDD_UP genes. The larger the font size, the greater the number of interacting partners per node. Fibronectin 1 (FN1), COL2A1 and CTNNB1 had the most interacting partners in this network.

Figure 2

Protein-protein interactions involving the three most connected proteins. Fibronectin 1 (FN1), COL2A1 and CTNNB1 were the first, second and third most connected proteins, with nine, eight and seven edges, respectively.

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 Thompson-based classification and the unsupervised molecular classification (p = 0.04799 for the 18,818-gene case and 0.09133 for the 1,984-gene case).

Figure 3

Unsupervised hierarchical clustering of 22 disc tissues of GSE23130. The pattern shown on the left was obtained by using all 18,818 genes while the result obtained by using only 1,984 genes with a minimum sample standard deviation of 1 is shown on the right.

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 up-regulated in punctured mouse tail discs that showed progressive degenerative changes and were fibrocartilaginous (Yang ). Earlier studies reported that the degradation fragments of fibronectin (Fn-f) accumulate in the disc during degeneration (Oegema ) and induce IVD degeneration in rabbits (Greg Anderson ). 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 ). The over-expression of these two genes was consistent with this hypothesis. Interestingly, fibronectin mRNA and protein levels were significantly down-regulated following the up-regulation of HTRA1 in rhesus monkey choroid-retina endothelial cells (RF/6A) and human umbilical vein endothelial cells (HUVECs) (Jiang ), 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 ). The relationship between COL2A1 gene polymorphisms and knee osteoarthritis were also investigated in Han Chinese women (Xu ).

Mutations in the COL11A1 gene (also known as STL2) have been identified in some people with Stickler syndrome (Martin ). Other mutations in this gene cause segments of DNA to be skipped when the protein is being made, resulting in an abnormally short pro-alpha 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 ). COL1A2 is a protein found in most connective tissues. Mutations in this gene are associated with asteo-genesis imperfecta, Ehlers-Danlos syndrome, idiopathic osteoporosis and atypical Marfan syndrome (Vasan ; Ward ; Hartikka ). However, the symptoms associated with mutations in this gene tend to be less severe than with mutations in the gene for α1 type I collagen since the α2 form is less abundant (Bou-Gharios ). Multiple messages for this gene result from multiple polyadenylation signals, a feature shared with most of the other collagen genes (Zuo ).

COL3A1 is a protein that, in humans, is encoded by the COL3A1 gene located on chromosome 2 (Cutting ). 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 ). 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 ), 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 ). 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 ). 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 ). Our results may provide clues for studying the biological mechanism of β-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.