An organ culture system to model early degenerative changes of the intervertebral disc II: profiling global gene expression changes

Introduction Despite many advances in our understanding of the molecular basis of disc degeneration, there remains a paucity of preclinical models which can be used to study the biochemical and molecular events that drive disc degeneration, and the effects of potential therapeutic interventions. The goal of this study is to characterize global gene expression changes in a disc organ culture system that mimics early nontraumatic disc degeneration. Methods To mimic a degenerative insult, rat intervertebral discs were cultured in the presence of TNF-α, IL-1β and serum-limiting conditions. Gene expression analysis was performed using a microarray to identify differential gene expression between experimental and control groups. Differential pattern of gene expression was confirmed using quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) or Western blot. Results Treatment resulted in significant changes in expression of more than 1,000 genes affecting many aspects of cell function including cellular movement, the cell cycle, cellular development, and cell death and proliferation. Many of the most highly upregulated and downregulated genes have known functions in disc degeneration and extracellular matrix hemostasis. Construction of gene networks based on known cellular pathways and expression data from our analysis demonstrated that the network associated with cell death, cell cycle regulation and DNA replication and repair was most heavily affected in this model of disc degeneration. Conclusions This rat organ culture model uses cytokine exposure to induce wide gene expression changes with the most affected genes having known reported functions in disc degeneration. We propose that this model is a valuable tool to study the etiology of disc degeneration and evaluate potential therapeutic treatments.


Introduction
Amongst the many potential pain generators in the lumbar spine, symptomatic disc degeneration is thought to be a significant contributor to low back pain (LBP) [1,2] and accounts for more than 25% of lumbar fusion surgery performed in the USA [3]. By some estimates, LBP results in direct and indirect health-care expenditures exceeding $100 billion per year, a considerable portion of which is due to discogenic pain [4]. Despite the scale of this clinical problem, many aspects of the pathogenesis associated with LBP remain incompletely characterized. The highly variable nature of associated symptoms and the presence of multiple potential sources of pain within the spine often confound efforts to accurately identify and study disc degeneration in live subjects. However, most authorities agree that there is a strong link between intervertebral disc degeneration and the LBP [5,6].
Establishing clinically relevant models of disc degeneration has also proven difficult. In contrast to humans, disc degeneration is not commonly seen in most quadrupedal animals, likely due to a variety of factors such as mineral composition, vertebral range of motion, extracellular matrix (ECM) composition, cell type, and weight distribution [7]. Different strategies have been developed to induce disc degeneration in animals including gene silencing, application of supraphysiologic loading, and disc injury. One theoretical advantage of using in vitro models to study disc degeneration is the ability to carefully control confounding environmental variables, which may differ between individual animals such as nutrition and the loading environment. To date, most widely used in vitro studies are performed with cells isolated from the intervertebral disc, which are then either grown in monolayer or suspended in a matrix that allows cells to assume a three-dimensional phenotype. To capture the benefits associated with in vitro models of disc degeneration while maintaining essentially normal cell density and distribution and the native ECM, we have developed and performed initial characterization of a disc organ culture model using intact rat intervertebral discs [8]. Since the hallmark catabolic processes such as matrix breakdown and decreased cellular biosynthesis during human degeneration are mediated by a number of cytokines, including IL-1β and TNF-α [9][10][11][12][13][14], we chose to mimic the degenerative state by treating organ-cultured discs with a cocktail containing both these cytokines.
Our initial investigations focused on targeted evaluation of a few key molecules known to be intermediaries in the degenerative cascade or essential components of the disc ECM [8]. In the present study, we hope to expand on our initial studies to demonstrate that this model replicates the degenerative phenotype in humans through analysis of global gene expression changes. Our results clearly indicate that this model is a good alternative to in vivo and cell-based in vitro models of intervertebral disc degeneration.

Disc harvesting and organ culture
The rat tissue samples were obtained under a protocol approved by the Institutional Animal Care and Use Committee (IACUC) of the Thomas Jefferson University (Protocol number 703H). Rat lumbar disc specimens were isolated using a method reported earlier by Ponnappan et al. [8]. Briefly, using a sterile surgical blade (number 15) the whole lumbar intervertebral discs (IVDs with the endplates were dissected (n = 6 discs/ animal) and maintained in organ culture in 12-well cell culture plates. All experimental lumbar discs were cultured in DMEM containing 1% FBS, 10 ng/ml IL-1β (R&D Systems, Inc., Minneapolis, MN, USA), 100 ng/ml TNF-α, 50 μg/ml L-ascorbate, 40 mM NaCl, antibiotics and antimycotics (Cellgro, Mediatech, Inc., Herndon, VA, USA). Control discs were cultured in DMEM containing 10% FBS, 50 μg/ml L-ascorbate, 40 mM NaCl and antibiotics without any cytokines. The discs were maintained in culture for 10 days. The complete medium was replaced every two days for both groups. In each experiment, six lumbar discs from one animal were used per group. The number of animals used for each study is indicated where appropriate.

Microarray analysis
Total RNA was extracted using the RNeasy Micro kit (Qiagen, Valencia, CA, USA, ) according to the manufacturer's instructions. DNase-treated RNA was quantified on a NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE, USA), followed by RNA quality assessment by analysis on an Agilent 2100 bioanalyzer (Agilent, Palo Alto, CA, USA). RNA amplification and labeling was performed from 50 ng total RNA by the WT-Ovation Pico RNA amplification system (NuGen Technologies, Inc., San Carlos, CA, USA) as described previously [15]. Each Affymetrix Rat Gene array (RaGene 1.0 ST v1) (Affymetrix, Santa Clara, CA, USA) was hybridized with fragmented and biotinlabeled target for 18 h. Arrays were then washed and stained using Gene-chip Fluidic Station 450, and hybridization signals were amplified using antibody amplification with goat IgG (Sigma-Aldrich, St Louis, MO, USA) and anti-streptavidin biotinylated antibody (Vector Laboratories, Burlingame, CA, USA). Chips were scanned on an Affymetrix Gene Chip Scanner 3000, using Command Console software. Background correction and normalization were done using Robust Multichip Average (RMA) with Gene-spring V 11.5 software (Agilent). A volcano plot was used to identify differentially expressed genes using the unpaired two-sample Student's t test, assuming significance for P values ≤ 0.05, and fold changes ≥1.5. A heat map was generated by Gene-spring software V 11.5. The differentially expressed gene list was loaded into Ingenuity Pathway Analysis (IPA-8.8) software to perform biological network and functional analyses following a previously described method [16]. In short, IPA was used to convert a list of differentially expressed genes into a set of relevant networks based on a continually updated database of known molecular pathways. This analysis is intended to suggest which cellular processes are most likely to be affected by the gene expression changes taking place in a given experiment. This process generates a network score that reflects the probability that a similar group of genes equal to or greater than the number in the network could be achieved through chance alone. Previous similar analyses have used a score of more than 10 as a cutoff for identifying relevant gene networks [16] with higher values corresponding to networks more likely to be affected by the expression changes seen in the experiment. Four animals for the control and four animals for the experimental group (24 discs/group) were used for this study, n = 4 independent experiments. The microarray data can be accessed at the Gene Expression Omnibus number GSE42611 [17].

RNA extraction and quantitative real-time PCR
The discs were collected and the annulus fibrosus (AF) was separated from the nucleus pulposus (NP) as in previously described method [8]. Total RNA was extracted from NP tissue (RNeasy Micro kit, Qiagen) according to the manufacturer's instructions. The RNA samples were treated with DNase I digestion prior to conversion into cDNA. Single-stranded cDNA was synthesized from 0.5 μg of total RNA using RNA to cDNA EcoDry Premix (Oligo dT) (Clontech Laboratories, Inc., Mountain View, CA, USA) and quantitative (q) RT-PCR was performed using SYBR Green PCR Master Mix (Applied Biosystems, Foster City, CA, USA) on an ABI 7900 HT sequence detection system. A melting curve analysis was performed to ensure primer specificity. The amount of PCR product was estimated using a relative standard curve quantification method. Expression was normalized by using housekeeping gene hypoxanthine phosphoribosyltransferase 1 (Hprt1). All the primers used were synthesized by Integrated DNA Technologies, Inc. (Coralville, IA, USA). Sequences for primers are presented in Table 1.

Senescence-associated b-galactosidase staining
The primary rat NP cells were isolated using a method reported earlier by Wang et al. [10]. Cells were plated in Lab-Tek II Chamber Slide System (Nalge Nunc International, Rochester, NY, USA) and treated with 10 ng/ml IL-1β and 10 ng/ml TNF-α for 10 days. The slides were rinsed with PBS, fixed with fixative solution for 15 min at room temperature, and incubated overnight at 37°C with fresh SA-β-gal staining solution (Senescence Detection Kit, Abcam). Cells were observed under a light microscope for development of blue color. The images were captured using a Nikon Eclipse E600 (Nikon, Tokyo, Japan).

Data analysis
Statistical analysis was performed using Sigma Plot 11.2 statistical software (Systat Software, Inc., Chicago, IL, USA). Differences between groups in continuous variables were analyzed by the Student's t test assuming significance at P < 0.05.

Microarray experiment
Following organ culture experiment ( Figure 1A), RNA was extracted and microarray analysis was performed. Figure 1B depicts the heat map of the significant NP gene expression profile changes between the control (Ctr) and experimental (Exp) groups, demonstrating  Volcano plot of microarray data depicting the significant differences (P < 0.05) in expression patterns between the control (Ctr.) and experimental (Exp.) groups, where the red plots represent a 1.5-fold or more difference in expression patterns between the two groups. (D) Pie chart showing differentially expressed genes based on major molecular and cellular functions. The total numbers of genes associated with each specific function are listed. Note that some genes may have multiple functions and can be classified in several categories (see Table 2), N = 4 independent experiments. more than 1,000 individual genes that were either upregulated or downregulated by more than 1.5-fold in the experimental model. In total, 1,036 genes demonstrated a significant upregulation or downregulation based on the microarray results. Figure 1C more clearly depicts the relative balance between genes either upregulated or downregulated by more than 1.5-fold (P < 0.05) in the experimental model compared to the control group using a volcano plot; significantly upregulated genes are depicted by red dots to the right of the center while significantly downregulated genes are denoted by red dots to the left of the center. From a functional standpoint, genes changing expression patterns in the experimental model could be divided into several different known cellular functions as seen in Figure 1D: 138 genes affect cellular movement, 73 genes affect the cell cycle, 156 genes affect cellular development, 139 genes are involved in cellular death and 180 genes are known intermediaries in cellular growth and proliferation. Genes contributing to these major cellular functions are shown in Table 2.
The 10 genes most highly upregulated and downregulated in the experimental group are depicted in Figure 2A and 2B, respectively, and the related expression levels are reported in Table 3. Upregulated genes with known roles in disc degeneration included interleukin-6 (Il-6), matrix metallopeptidase 3 (Mmp3), and a disintegrin and metalloproteinase with thrombospondin motifs 5 (Adamts5). Highly downregulated genes with known roles in maintenance of normal disc phenotype or as markers for NP cells included Krt19, asporin (Aspn), collagen I alpha-1 (Col1a1), collagen I alpha-2 (Col1a2) and insulin-like growth factor 1 (Igf1).

RNA extraction and quantitative real-time PCR Validation of differentially expressed genes by real-time RT-PCR
Confirmation of microarray data with qRT-PCR was performed for selected genes from the top 10 upregulated and downregulated genes as shown in Table 3. Also included in these confirmatory experiments was Sdc-4, as it was also significantly upregulated in the array data and has recently been shown to be a key regulator of Adamts5 activity in the IVD [10]. Experimental group demonstrated significant upregulation of chemokine (C-C motif) ligand 20 (Ccl20), Il-6, Adamts5, and Mmp3 by almost 216-fold (P < 0.001), 45-fold (P = 0.003), 26-fold (P = 0.02) and 42-fold (P = 0.048), respectively, compared to control ( Figure 2C and 2D). Alternatively, Krt19 and Krt8 were significantly downregulated in the experimental group by 27-fold (P < 0.001) and 17-fold (P < 0.001), respectively, compared to controls ( Figure 2E). Col1a1 and Col1a2 were not evaluated using qRT-PCR in this study as these genes have previously been shown to be downregulated in experimental discs using the same organ culture model [8].
Notably, the relative changes in gene expression between experimental and control groups were larger when analyzed by qRT-PCR suggesting that microarray analysis may be a conservative evaluation of the gene expression differences between the control and the experimental discs.

Western blot analysis
Western blot analysis was performed to demonstrate the effect of treatment on the expression of SDC-4 and KRT19. Similar to the results for the microarray and the qRT-PCR experiments, Western blot analysis showed a robust downregulation of KRT19 protein in the experimental group compared to the control group ( Figure  2F). There was a strong increase in SDC4 levels in the experimental discs when compared to the control discs ( Figure 2F).  Figure 3A is associated with cell death, cell cycle, and DNA replication, recombination and repair, and has a score of 31. The network depicted in Figure 3B is associated with cellular growth and proliferation, tissue development and cancer, and has a score of 27. The network shown in Figure 4A is associated with cellular development and vascular function and development, and has a score of 27. The network illustrated in Figure 4B is associated with carbohydrate metabolism and protein synthesis and has a score of 16. Finally, the network in Figure 5 depicts genes related to cellular movement, immune cell trafficking, cell-to-cell signaling and cellular interactions and has a score of 24.

Effect of cytokines on NP cell senescence
Staining for SA-β-gal was used to measure the effect of proinflammatory cytokines on phenotype of primary rat NP cells. Increased blue staining was evident in the cytoplasm of the cells treated for 10 days compared to the control group, indicating onset of senescent cell phenotype ( Figure S1 in Additional file 1).
The phenotypic and genotypic changes associated with disc degeneration are complex, reflecting the interrelatedness of many distinct processes. We found changes in the expression of genes involved in a number of cellular functions; genes controlling cellular growth and proliferation (181 genes) and cellular development (156 genes) were the most frequently affected. Similarly, IPA network analysis demonstrated the highest gene network score amongst genes involved in cell death, cell cycle, and DNA replication, recombination and repair. This finding, and the increased incidence of SA-β-gal in cytokine-treated NP cells, suggests that during degeneration disc cells may be particularly susceptible to processes such as senescence [28], at the same time as identifying genes that are likely to be related to these processes. Furthermore, these IPA networks highlight the web of interrelated pathways that are differentially regulated in discs undergoing degenerative changes; not only may changes in expression of a single gene have wide downstream consequences but these figures hint that the existence of a single key mediator of disc degeneration is unlikely.
Analysis of the top 20 most affected genes indicate that many genes are known mediators of disc degeneration or serve a role in maintaining normal disc function. Prior characterization of gene expression patterns from degenerative human IVD tissue has shown a similar breadth of changes in global expression patterns [29][30][31]. These gene expression analyses suggest that many different processes are altered in disc degeneration, although most current research tends to focus on the patterns of inflammatory mediator expression. In agreement with recent studies that showed enhanced chemokine-mediated macrophage migration by inflammatory cytokines [32], one of the highly upregulated molecules identified in this study is Ccl20. This chemokine is a major regulator of macrophage and mast cells migration and has been reported to be upregulated in other arthritic disorders [33,34]. Gruber et al. compared global gene expression profiles between discs with different degrees of degeneration and found broad changes in gene expression patterns, both with respect to the large number of affected genes and the variety of known functions [29]. In contrast to the present study, few of the 47 genes identified by Gruber et al. to be differentially expressed had known roles in disc degeneration. In another study evaluating differences between disc tissue harvested immediately post-mortem from patients without known spinal pathology and tissue harvested from patients with degenerative discs, Zhang et al. [30] reported more than 500 differentially expressed genes. Consistent with the results from our study, Zhang and colleagues reported a significant upregulation of two known regulators of disc degeneration: Mmp3 (3-fold) and Il-6 (2.5-fold). Finally, Gruber et al. [31] recently utilized microarray analysis and found more than 4,500 differentially expressed genes between early and advanced stages of disc degeneration. Moreover, the authors identified 23 with known roles in the homeostasis of the ECM including collagen subunits, Aggrecan, Adamts and metalloproteinase inhibitor 3 (Timp3).
Certain genes identified in our screen have previously been studied for their role in the pathobiology of intervertebral disc degeneration. Moreover, the similarity with the previously identified expression patterns in human disc degeneration suggests that our model accurately recapitulates natural disc degeneration. IL-6 is an inflammatory cytokine that has been associated with disc herniation and discogenic pain in clinical studies [35,36] but has been less well characterized with respect to its cellular function compared to other inflammatory mediators such as IL-1 or TNF-α. Studer et al. [37], recently described contribution of IL-6 in the presence of soluble IL-6 receptor (sIL-6R) in amplification of the inflammatory response to IL-1 and TNF-α in NP cells. The authors suggested that in the NP, IL-6/sIL-R pathway may play a role in potentiating specific functions of IL-1 and TNF-α such as downregulation of matrix production and upregulation of prostaglandins and MMP3. We found that both Mmp3 and Adamts5 were highly upregulated in our model. MMP3 is an extracellular zinc-dependent proteinase involved in digestion of noncollagen matrix proteins and regulates ECM homeostasis in healthy discs. Upregulation of Mmp3, however, is associated with excessive degradation of the ECM, ultimately depleting the hydrophilic nature of the NP and the associated ability to distribute compressive load. Similarly, the ADAMTS peptidases that degrade important proteoglycans including aggrecan [10,38] are shown to be upregulated in degenerative discs by both Le Maitre et al. [39] and Pockert et al. [40]. SDC4 is a cell surface proteoglycan that has been shown to be essential in activation of the ADAMTS5 within the intervertebral disc, thereby serving as an important regulator of ECM homeostasis [10]. Our study found substantial upregulation of Sdc4, confirming the important role Adamts play in ECM regulation during disc degeneration. It is important to comment that there is a strong correlation between expression of senescent biomarkers and increased gene expression of catabolic-degrading enzymes in the disc. Thus up-regulation of many of the catabolic enzymes in our model is likely a reflection of senescent cell phenotype as indicated by our SA-β-gal staining [26].
Disc degeneration also affects synthesis of ECM proteins. Collagen I is an important structural protein within the disc, particularly in the AF, which must resist high tensile loads to maintain the disc shape and spine alignment during axial loading. Disc degeneration is typically associated with an upregulation of collagen I that leads to a loss of compliance and hardening of the NP. Interestingly, we found downregulation of two collagen I alpha chains, contrary to what would be expected. Although downregulation of collagen II and upregulation of collagen I is eventually seen in a degenerative IVD, one explanation for these findings is that the initial reparative attempts by the IVD normalize ECM expression patterns, which only later become irrevocably altered.
Downregulation of asporin in the treated cells was also contrary to expectations. Asporin was found to be expressed at higher levels in human discs with more advanced degeneration at both the protein and mRNA level [41]. Furthermore, expression of a particular asporin allele (D14) has been associated with the development of disc degeneration in both Chinese and Japanese populations [42]. Little, however, is known about the molecular basis for these findings as the role of asporin in the IVD has not been thoroughly characterized. Further understanding of the role of this molecule in disc degeneration may provide greater insight into potential explanations for downregulation of asporin in organ culture models.
Other researchers have described organ culture models to study disc degeneration, but used different processes to induce degeneration. Roberts et al. [18], Chen et al. [19] and Jim et al. [43] have all recently described systems that utilize injections of chemolytic enzymes to degrade the ECM and induce a degenerative phenotype. Although these models result in a phenotype consistent with disc degeneration, the homology with the human intervertebral disc degenerative process is questionable as there is little evidence that the nonspecific global chemolysis is the driving force behind human IVD degeneration [29,30,43]. More work is needed to establish whether the use of such enzymes is an appropriate way to induce structural degeneration and inflammation for use in modeling human disc degeneration. Similarly, many organ culture models are based on needle puncture and may be limited by this mechanism of disc injury. Annular needle puncture in animal models and humans results in clinically significant iatrogenic disc degeneration, however, it is unclear whether this model is an accurate representation of disc degeneration commonly seen in humans [44,45] that does not occur via this mechanism. In contrast, induction of disc degeneration using IL-1β and TNF-α mimics the pathways that drive inflammatory processes and disease phenotype in humans.

Conclusions
Our results clearly demonstrate a widespread differential regulation of genes largely consistent with previous research into the molecular basis of disc degeneration. Furthermore, results from our study found differential expression in many genes known to be mediators of inflammatory processes associated with disc degeneration and regulators of extracellular matrix production. Taken together, these findings suggest that the rat organ culture model described can be used for further studies designed to better understand the etiology of disc degeneration and evaluate potential therapeutic treatments.

Additional material
Additional file 1: Figure S1. Senescence-associated β-galactosidase staining of rat nucleus pulposus (NP) cells following treatment with TNFα and IL-1β for 10 days. The result indicates that the number of SA-βgal-positive NP cells is increased in the experimental (B) versus the control group (A). Authors' contributions DZM and SA carried out the experimental work, analyzed the data and drafted the manuscript. CKK helped design the studies, secured funding and co-wrote the manuscript. HBM, ARV, IMS, DGA, and TJA helped design the study and prepare the final manuscript. MVR designed the study, helped analyze the data, secured funding and co-wrote the manuscript. All authors read and approved the final manuscript.