Transcriptional profiling of intervertebral disc in a post‐traumatic early degeneration organ culture model

Abstract Introduction The goal of this study is to characterize transcriptome changes and gene regulation networks in an organ culture system that mimics early post‐traumatic intervertebral disc (IVD) degeneration. Methods To mimic a traumatic insult, bovine caudal IVDs underwent one strike loading. The control group was cultured under physiological loading. At 24 hours after one strike or physiological loading, RNA was extracted from nucleus pulposus (NP) and annulus fibrosus (AF) tissue. High throughput next generation RNA sequencing was performed to identify differentially expressed genes (DEGs) between the one strike loading group and the control group. Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes analyses were performed to analyze DEGs and pathways. Protein‐protein interaction (PPI) network was analyzed with cytoscape software. DEGs were verified using qRT‐PCR. Degenerated human IVD tissue was collected for immunofluorescence staining to verify the expression of DEGs in human disc tissue. Results One strike loading resulted in significant gene expression changes compared with physiological loading. In total 253 DEGs were found in NP tissue and 208 DEGs in AF tissue. Many of the highly dysregulated genes have known functions in disc degeneration and extracellular matrix (ECM) homeostasis. ACTB, ACTG, PFN1, MYL12B in NP tissue and FGF1, SPP1 in AF tissue were verified by qRT‐PCR and immunofluorescence imaging. The identified DEGs were involved in focal adhesion, ECM‐receptor interaction, PI3K‐AKT, and cytokine‐cytokine receptor interaction pathways. Three clusters of PPI networks were identified. GO enrichment revealed that these DEGs were mainly involved in inflammatory response, the ECM and growth factor signaling and protein folding biological process. Conclusion Our study revealed different DEGs, pathways, biological process and PPI networks involved in post‐traumatic IVD degeneration. These findings will advance the understanding of the pathogenesis of IVD degeneration, and help to identify novel biomarkers for the disease diagnosis.


| INTRODUCTION
Low back pain (LBP) is a significant healthcare concern worldwide, which leads to major economic and social burdens that affect millions of people and is a leading reason for seeking medical attention with the lifetime prevalence approaching 84%. 1 Degeneration of the intervertebral disc (IVD) is a major contributing factor to chronic LBP. The presence of individuals over the age of 50 with evidence of intervertebral disc degeneration (IVDD) approaches 90%, correlating with lifetime prevalence of LBP. 2 IVDD is characterized by increased imbalance between extracellular matrix (ECM) breakdown and matrix synthesis. 3,4 This imbalance results in dehydration of the central nucleus pulposus (NP), reduced proteoglycan content, decreased cellularity, diminished endplate density, and disruption of the annulus fibrosus (AF). 5 Environmental exposures, as well as genetic and epigenetic factors have been associated with disc degeneration and altered ECM synthesis in disc tissues. 6 Many complex and interdependent factors, including aging, decreased nutrient supply, unphysiological loading (e.g., overloading), have been implicated in the initiation and progression of IVDD. Among these factors, unphysiological mechanical loading is generally thought to play a crucial role in the process of IVDD, which is associated with ECM degradation and cellular loss. 7 Different methods have been developed to induce IVDD in animals including gene silencing, application of supraphysiological loading, and disc injury. The advantage of using ex vivo models to study IVDD is the ability to carefully control confounding environmental variables, which may differ between individual animals such as nutrition and the loading environment. 8 IVD whole organ culture with different degeneration initiating factors can mimic various phenotypes of disc degeneration. 9 The ability to treat symptomatic IVDD effectively is hindered by an incomplete understanding of the biological processes that control IVD development, function, and disease. Consequently, the clinical management of IVD pathologies remains very limited, with no options at present for early intervention or predictive patient screening.
Therefore, there is a need for an improved understanding of the molecular pathophysiological mechanisms underlying IVDD, which is essential for diagnosis and the development of novel therapeutic approaches. Recently, the molecular basis of IVDD has received increased attention in research, which has substantially improved the understanding of the biology underlying this process. Studies investigating the molecular changes associated with the pathophysiology of IVDD have established criteria to distinguish degenerative IVDs. 4,10 Studies evaluating transcriptome data using microarrays have provided us with an initial understanding of the molecular mechanisms underlying disc biology. 11,12 Furthermore, high-throughput screening of human patient samples may identify potential biomarkers of IVDD, leading to more precise diagnostic criteria, classification of disease progression, and prognosis.
Genes work in synergy with each other to perform biological functions. They simultaneously interact with multiple genes and trigger a variety of changes that lead to diverse reactions. The functional annotation of regulated genes, using Gene Ontology (GO), has enabled the identification of severely affected groups of genes that correlate with the disease phenotypes. 13 Therefore, the analysis of the gene expression profile by applying bioinformatics methods remains necessary to identify differentially expressed genes (DEGs) in IVDD and further elucidate the potential pathogenesis mechanisms of the disease.
High impact loading is one of the major causes leading to disc herniation. 14,15 The herniation may not occur right after the one strike overload, but years after. Therefore, the one strike overload model was selected to mimic the post-traumatic IVD degeneration. [16][17][18] Previously, we have established an IVDD model using one strike loading to mimic the post-traumatic pathological changes in whole organ cultured IVD. A single hyperphysiological mechanical compression applied to healthy bovine IVDs caused significant drop of cell viability, altered the mRNA expression in the IVD, and induced ECM degradation. 19 The present study aimed to identify the DEGs induced by one strike loading and further analyze their functions and pathways associated with the progression of IVDD by utilizing bioinformatics methods. To achieve this objective, we evaluated RNA profiles of NP and AF cells from the bovine IVD organ culture model utilizing high throughput next generation RNA sequencing. Furthermore, immunofluorescent staining was used to validate the expression of selected DEGs in human degenerated IVD tissue. Caudal bovine IVDs (n = 8) were obtained from 4 to 10 months old animals as previously described. 20,21 Briefly, after removal of the soft tissues, IVDs comprising endplates were harvested using a band-saw.
The endplates (1-2 mm thickness) were rinsed with Ringer solution using the Pulsavac jet-lavage system (Zimmer, Warsaw, Indiana). The discs were further incubated in phosphate buffered saline (PBS) solution with 1000 units/mL penicillin and 1000 μg/mL streptomycin for 10 min. The cleaned discs were transferred to 6 well-plates and cultured in an incubator at 37 C, 85% humidity and 5% CO 2 until the next day. The culture medium was composed of Dulbecco's modified Eagle medium (DMEM) containing 4.5 g/L glucose and supplied with 10% fetal calf serum, 100 units/mL penicillin, 100 μg/mL streptomycin (all products from Gibco, Basel, Switzerland) and 0.1% Primocin (Invitrogen, San Diego, California). Discs (mean IVD height 9.47 ± 1.27 mm, mean IVD diameter 15.96 ± 1.75 mm) were randomly assigned to one of two groups: one strike loading model (n = 8/ group), physiological loading control (n = 8/group).

| Loading protocol
After dissection on day 0, the IVDs were cultured under physiological loading at 0.02-0.2 MPa, 0.2 Hz, 2 hours per day within a custom designed bioreactor on day 1. 9,21,22 In our previous study, we found that this loading protocol caused a disc height loss around 10%, 21 which is similar to the physiological disc height loss in human lumbar discs after diurnal activities. 23 The bioreactor was maintained in an incubator at 37 C, 85% humidity and 5% CO 2 . On day 2, IVDs from the one strike loading group were subjected to one strike loading in a custom designed chamber or physiological loading within the bioreactor; the physiological loading control group received physiological loading. One strike loading was applied using a Mini Bionix 858 MTS machine. Isolated discs were held under a mild load of 10N for 3 min to reach steady contact with the load cell. Discs were then subjected to a single compression at 50% of disc height at a velocity of 10 mm/s. The average peak stress achieved in the one strike loading model group was 35.461 ± 8.274 MPa. The culture medium was composed of DMEM

| RNA extraction
Disc tissue was collected on day 3 for gene expression measurement.
Cartilaginous endplates of each IVD were removed, and NP was harvested using a biopsy punch and the outer AF tissue was harvested from the outer region of AF using a scalpel blade.
For RNA extraction, approximately 150 mg of each NP and AF tissue was used. Tissue samples were digested in 2 mg/mL pronase for 1 h at 37 C, flash frozen, pulverized in liquid nitrogen, and homogenized using a TissueLyser (Qiagen, Venlo, Netherlands). 24 Total RNA was extracted with TRI Reagent (Molecular Research Center).

| RNA sequencing
In total 12 samples were performed with RNA sequencing, which included 3 NP samples and 3 AF samples from the control group, and 3 NP samples and 3 AF samples from the one strike model group. Library preparation and RNAseq were performed at the Genomics Core Facility In brief, 10 ng of total RNA was used to generate amplified cDNA with priming at the 3 0 end as well as randomly throughout the transcriptome. Next, libraries were prepared by fragmentation of the ds-cDNA, end repair to generate blunt ends, and adaptor ligation. Finally, the rRNA AnyDeplete probe mix was added to facilitate degradation of unwanted ribosomal RNA transcripts. The remaining library molecules were then PCR amplified (6 cycles) to produce   Table 1. Comparative Ct method was performed for relative quantification of target mRNA with Ribosomal protein large P0 (RPLP0) as endogenous control. 20

| Preprocessing and differential analysis
Raw data from RNA sequencing were converted into a recognizable format with the package affy of R (http://bioconductor.org/packages/ release/bioc/html/affy.html,version 1.54.0), and missing values were then inferred by a method based on k-nearest neighbors (k-NN). Following background correction and data normalization with the median method, differential analysis between one strike loading samples and controls was performed using limma package (version 3.32.5). For statistical analysis and the assessment of differential expression, an empirical Bayes method from limma was employed to moderate the standard errors of the estimated log2-fold changes. The basic statistic used for significance analysis is the moderated t-statistic, which is computed for each probe and for each contrast. Moderated t-statistics lead to P-values in the same way as ordinary t-statistics, except that the degrees of freedom are increased, reflecting the great reliability associated with the smoothed standard errors. Limma functions top table and decide tests were used to summarize the results of the linear model, perform hypothesis tests and adjust the P-values for multiple testing. The results obtained include log2-fold changes, standard errors, t-statistics, and P-values. 25 Logj2(fold change)j >1 and P <0.05 were set as the cut-offs to screen out DEGs.

| GO functional enrichment analysis of DEGs
In order to understand the importance of genes and to identify dis-

| Pathway enrichment analysis of DEGs
The Kyoto encyclopedia of Genes and Genomes (KEGG: http://www. genome.jp/kegg) database is a collection of online databases consisting of genomes, enzymatic pathways, and biological chemicals. In our study, KEGG pathway enrichment analysis was performed to determine the function of DEGs using KOBAS 2.0 with a threshold of P <.05. KOBAS 2.0 is a web server that provides a comprehensive functional annotation tool. The biological pathways were associated with genes based on mapping to genes with established annotations.
Statistical testing was performed to identify statistically significantly enriched pathways and diseases. 27   network. Significant modules were identified according to the clustering score using the following criteria: "Degree cutoff = 2," "node score cutoff = 0.2," "Haircut = true," "Fluff = false," "k core = 2," and "max depth = 100." The clustering modules having high node scores and connectivity degrees were considered as biologically significant clusters.

| PPI network modules
To reveal associations among genes within the most enriched network module, the Cytoscape Web plugin was used to visualize the input genes in a network graph. STRING (http://string-db.org) was used to analyze the gene related protein interaction network. The names of genes related to PPI network are shown in Figure 5. KRT8, KRT18, CDH2, and CCND1 were identified as a subset of hub genes in the PPI network in NP tissue, which were related to aging and proliferation. IL6, CXCL6, CCR7, and CCL20 were identified as a subset of hub genes in AF tissue which play an important role in the inflammatory response of the IVD. COL2A1, COL9A2, COL11A2, COL19A1, and COL27A1 were identified as a subset of hub genes in AF tissue related to ECM composition.

| Validation of DEGs by real-time qRT-PCR
Confirmation of gene sequencing data with qRT-PCR was performed for selected genes from the DEGs. These genes are related to ECM- annulus fibrous tissue of degenerative and nondegenerative discs. 31 The result revealed that 238 DEGs was dysregulated in various cellular functions including cell proliferation and inflammatory response. 31 Guo et al 32 used microarray to compare the gene expression profiles between patients with IVD degeneration and control. A total 93 and 114 DEGs were identified in NP and AF tissue, respectively. GO analysis showed that the DEGs may be involved in various processes, including cell adhesion, biological adhesion, and ECM organization. 32 KEGG analysis showed that the DEGs was involved in focal adhesion and the p53 signaling pathway. 32 Wang et al 33  SPP1 could activate the integrin proteins in NP cells and has close connections with COL11A1, COL3A1, and COL2A1. 37 Since integrins are a class of cell adhesion molecules that regulate interactions between cells and their surrounding matrix, activation of integrin alpha (ITGA) will assist NP cells to interact with the ECM, specifically collagen molecules and then to potentially restore the IVD function. 38 Marfia et al reported greater expression of SPP1 in degenerative IVD compared to herniated IVD, and SPP1 was only detected in degenerative IVD tissue. It is speculated that SPP1/OPN may be a marker for the severity of disc degeneration. 39 The actin cytoskeleton is essential for the proper functioning of many cellular processes, including maintenance of the cell shape, chemotaxis, cell movement, adhesion, transport of cellular organelles, mitosis, replication, transcription, and even DNA repair. 40,41 Studies had identified that β-actin predominates in stress fibers while γ-actin is enriched at the leading edge of cell surface. 42 Studies with immature bovine IVDs showed that AF cells respond to cyclic tensile strain (10%, 1 Hz) with increase in the expression of β-actin and β-tubulin at both the transcriptional and protein level. 43 Cell-matrix adhesion has essential roles in a number of important biological processes, including cell motility, proliferation and differentiation, and the regulation of gene expression and cell survival; at contact points between the cell and ECM, specialized structures termed focal adhesions are formed. 44 The response of cells to loading, and the specific mechanotransduction pathways mediating responses to load, are dependent on cell morphology, cell-cell interactions, and cell-ECM interactions. Loading of cells can act directly on cytoskeleton networks by promoting F-actin reorganization, which can act as a mechanism for transducing mechanical stimuli to the nucleus. 43,45 Actin cytoskeletal remodeling is largely mediated by the RhoA/ROCK signaling pathway, which is one important part of focal adhesion pathway also identified by our research. The RhoA/ROCK signaling pathway is involved in many responses such as cell migration, polarization, stress fiber, and adhesion formation. 46 RhoA/ROCK signaling is also altered in response to mechanical stimuli, which creates complex interactions between morphological and biomechanical changes. 47 Supporting these findings, Guo et al 32  An inflammatory response is thought to initiate IVDD, and proinflammatory molecules, secreted by IVD cells are considered to mediate IVDD. 48 These cytokines trigger a range of pathogenic responses by the disc cells that can promote autophagy, senescence and apoptosis. 6,49,50 The resulting imbalance between catabolic and anabolic responses leads to degeneration, as well as herniation and radicular pain. 48 In the present study, IL6, CXCL6, CCR7, and CCL20 were identified as a subset of hub genes in the PPI network which was constructed in the AF tissue. Accordingly, it may be speculated that this combination of genes has a key role in the immediate response of AF cells to mechanical stress. IL6 plays a critical role in IVDD. 51  released in the culture medium from the one strike model was significantly higher than the control group. 19 The results showed the disturbance of the balance between anabolism and catabolism of the ECM during the IVD degeneration. 19 Disc degeneration also affects the synthesis of ECM proteins. Collagen and aggrecan are important structural proteins and proteoglycans within the disc, which must resist high tensile loads to maintain the disc shape. These ECM protein genes are often dysregulated during IVDD. For example, disc degeneration is typically associated with an upregulation of collagen I that leads to a loss of compliance and hardening of the NP. Downregulation of collagen II and upregulation of collagen I in NP is eventually seen in a degenerative IVD. One explanation for these findings is that the initial reparative attempts by the IVD restore ECM expression patterns, which only later become irrevocably altered. 8  Hwang et al 68  There are several limitations in our study. First, the main focus was to investigate the mechanism of early degeneration change induced by mechanical loading. However, IVDD can also be induced by other factors such as nutrient insufficiency, cell senescence, and these factors may have a significant influence on the biomolecular signaling events within IVD degeneration. Second, follow up of the subsequent degeneration mechanism was not investigated in the current study. Third, immunofluorescence staining was performed on degenerative human IVD samples to verify the protein expression of selected DEGs. However, this is only preliminary investigation without healthy human IVD tissue as control group. Further research on the expression level of the DEGs in healthy control and human IVD tissue at different degeneration grades will be needed to verify their potential as biomarker or therapeutic target for post-traumatic early IVD degeneration. Finally, only a set of genes that were already known to be related with IVD degeneration were validated by qRT-PCR and immunofluorescence staining.

| CONCLUSIONS
In conclusion, the present study investigated the early transcriptome changes after single high impact load of the IVD via bioinformatics analysis. A total of 252 DEGs were identified in the impacted NP and 208 DEGs were identified in the impacted AF. The combination of ACTB, ACTG, PFN1, MYL12B, and SPP1 may serve as biomarkers for early responses to mechanical overloading. These DEGs may be involved in focal adhesion, ECM-receptor interaction, PI3K-AKT and cytokine-cytokine receptor interaction pathways during the early phase of IVD degeneration caused by mechanical loading. GO enrichment revealed that these DEGs were mainly involved in inflammatory response, the ECM, and growth factor signaling and protein folding processes.
These molecules are potential novel targets of therapeutic strategy for the treatment or prevention of post-traumatic IVD degeneration. The findings shed new light on the pathophysiology of IVD degeneration and may guide rational treatments, including gene therapy, cell therapy, and tissue engineering approaches for IVD degeneration and LBP.