Introduction
Intervertebral discs lie between adjacent vertebrae in the spine, forming a fibrocartilaginous joint to allow slight movement of the vertebrae (1). The layers of fibrocartilage contain the nucleus pulposus, which functions as a shock absorber dissipating compressive forces outward from its center to the surrounding annulus fibrosus (2). However, the anulus fibrosus becomes weaker with increasing age (3), resulting in intervertebral disc degeneration (IDD). Beyond age 40, >60% of individuals show symptoms of IDD; furthermore, IDD is the most common cause of disability among workers aged 18–64 years (4). IDD is characterized by decreases in intervertebral disc function and height due to cell loss through apoptosis, increased breakdown of matrix or altered matrix synthesis, with the underlying pathological processes being complex (5–7).
Numerous studies have intended to explore the molecular mechanisms involved in IDD. Vo et al (4) reported that IDD was a consequence of increased catabolism of the extracellular matrix (ECM), since the proteolytic degradation of ECM macromolecules led to marked structural changes of the intervertebral disc. These catabolic processes are mediated by a number of cytokines in the nucleus pulposus, among which interleukin (IL)-1β and tumor necrosis factor (TNF)-α have been suggested to have crucial roles in the development of IDD. However, the involvement of IL-1β and TNF-α in IDD has remained to be fully elucidated.
Markova et al (8) cultured rat intervertebral discs in the presence of IL-1β, TNF-α and serum-limiting conditions to mimic a degenerative insult to identify the differentially expressed genes (DEGs) between experimental and control groups. Their microarray data (no. GSE42611) have been deposited at the National Center of Biotechnology Information Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) and were used in the present study for a bioinformatics analysis. DEGs between experimental and control samples were identified, Gene Ontology (GO) and pathway enrichment analyses were performed and a protein-protein interaction (PPI) network was constructed. The present study aimed to identify the DEGs and pathways associated with IDD caused by the presence of IL-1β and TNF-α.
Materials and methods
Affymetrix microarray data
The microarray expression profile dataset GSE42611 (8) was downloaded from the GEO (http://www.ncbi.nlm.nih.gov/geo/) database based on the platform of Affymetrix Rat Gene 1.0 ST Array [transcript (gene) version; Affymetrix Inc., Santa Clara, CA, USA]. The dataset contained four experimental nucleus pulposus (ENP) samples and four control nucleus pulposus (CNP) samples. The experimental lumbar discs had been cultured in Dulbecco's modified Eagle's medium (DMEM; R&D Systems, Inc., Minneapolis, MN, USA) containing 10 ng IL-1β (R&D Systems, Inc.), 100 ng/ml TNF-α (R&D Systems, Inc.), 1% fetal bovine serum (FBS; R&D Systems, Inc.), 50 µg/ml L-ascorbate (Cellgro; Corning Incorporated, Corning, NY, USA), 40 mM NaCl (Cellgro; Corning Incorporated), antibiotics and antimycotics (Cellgro; Corning Incorporated), while the control discs were cultured in DMEM containing 10% FBS, 50 µg/ml L-ascorbate, 40 mM NaCl, antibiotics and antimycotics.
Data preprocessing and differential expression analysis
The original array data were converted into gene symbols and then subjected to background correction and quartile data normalization using the robust multiarray average (9) algorithm in the oligo package (10), which is available through BioConductor (http://www.bioconductor.org).
The paired Student's t-test based on the Limma package (11) in R was used to identify DEGs between ENP and CNP samples. Multiple testing correction was performed using the Benjamini-Hochberg method (12) to obtain the adjusted P-value. Subsequently, the log2-fold change (log2FC) was calculated. Only genes with an adjusted P<0.05 and a |log2FC|>1.0 were regarded as DEGs.
GO and pathway enrichment analyses
GO (http://www.geneontology.org/) (13) is a database used for unification of biological data, which comprises a structured, defined and controlled vocabulary for large-scale gene annotation. The Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.genome.jp/) (14) database is a collection of online databases of genomes, enzymatic pathways and biological chemicals.
The present study performed GO and KEGG pathway enrichment analyses to determine the function of DEGs using the Database for Annotation, Visualization and Integrated Discovery (http://david.abcc.ncifcrf.gov/) (15) online tool, which is a comprehensive functional annotation tool for associating functional terms with gene lists using a clustering algorithm. P<0.05 and a gene count >2 were set as thresholds.
PPI network construction
The Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org/) (16) database is a pre-computed global resource, which was designed to predict functional associations between proteins. In the present study, the STRING online tool was applied to analyze the PPIs of DEGs with experimentally validated interactions with a combined score of >0.4 considered significant.
The previously obtained biological networks indicated that the majority of the biological networks are characterized as scale-free (17). Thus, the connectivity degree was assessed by statistical analysis of networks to identify important nodes, which are being referred to as hub proteins (18).
Prediction of regulatory networks of transcription factors (TFs)
The TRANSFAC® database (19), which is available at http://www.gene-regulation.com, consists of TFs, their target genes and regulatory binding sites. In the present study, the transcriptional regulatory network within the PPI network was predicted based on the TRANSFAC® database.
Network module identification and functional enrichment analysis
After integrating the transcriptional regulatory network into the PPI network, the modular complex detection plugin (20) in cytoscape (http://www.cytoscape.org/) (21) was used to identify the network modules. The selected modules with high node scores and connectivity degrees were the subjected to GO and KEGG pathway enrichment analyses.
Results
Identification of DEGs
After data pre-processing, a total of 536 DEGs were obtained between ENP and CNP samples. Among these DEGs, 246 were upregulated and 290 were downregulated.
GO and pathway enrichment analyses
The top two clustering groups obtained by GO enrichment analysis are shown in Table I. The upregulated DEGs were mainly enriched in biological processes (BPs) associated with inflammatory response and responses to organic substances. The downregulated DEGs were enriched in the BP terms associated with cell adhesion and collagen fibril organization.
The results of the pathway enrichment analysis are shown in Table II. The upregulated DEGs were enriched in nine pathways, which included apoptotic pathways and the NOD-like receptor signaling pathway. The downregulated DEGs were also enriched in nine pathways, including pathways of ECM-receptor interaction and focal adhesion.
PPI network construction
Based on the STRING database, 1,345 PPI pairs were obtained and the PPI network was constructed (Fig. 1). The PPI network contained 36 DEGs with a degree of >25, referred to as hub genes, which included IL6 (degree, 67), COL1A1 (degree, 36), NFKB1 (degree, 31) and HIF1A (degree, 26).
Regulatory networks of TFs
Based on the PPI network and the TRANSFAC® database, 81 PPI pairs of transcriptional regulatory interactions were obtained. In addition, seven TFs were obtained, namely RAR-α, ANPEP, ETS2, ATF3, EGR1, HIF1A and NFKB1.
Network module identification and functional enrichment analyses
Using the MCODE plugin in cytoscape, four modules were obtained. Module 1 contained 29 nodes and 376 edges; however, it only contained downregulated DEGs and no TFs. Module 2 contained 20 nodes and 74 edges, and included the upregulated TFs HIF1A and NFKB1 (Fig. 2). Module 3 contained 10 nodes and 33 edges while module 4 contained 14 nodes and 42 edges. Module 2 was further subjected to GO and pathway enrichment analyses, and the results were shown in Table III. The DEGs in module 2 were predominantly enriched in GO terms associated with the regulation of cell proliferation and wound healing. In addition, they were also enriched in pathways associated with cancer and focal adhesion.
Discussion
IDD is present in adults with degenerative disc disease and is considered a major source of back pain in middle-aged adults, leading to a decrease of life quality or even disability (22). To date, the molecular mechanisms involved in the pathology of IDD have not been fully elucidated. In the present study, a total of 246 upregulated and 290 downregulated DEGs were identified between ENP and CNP samples. Several upregulated DEGs, including NFKB1, were enriched in apoptotic pathways, while downregulated DEGs, including COL1A1, were enriched in the pathway of ECM-receptor interaction. In addition, in the PPI network, IL6, COL1A1, NFKB1 and HIF1A were hub genes with a high connectivity degree. Furthermore, module analysis indicated that the DEGs in module 2 were predominantly enriched in GO terms associated with cell proliferation and pathways associated with cancer. The results suggested that these genes and pathways may be candidates for IDD diagnosis and treatment.
A central feature of IDD is loss of tissue cellularity, which may be the result of programmed cell death (23). The results of the pathway enrichment analysis indicated that apoptosis was the most significant pathway. Apoptosis is a type of programmed cell death and is characterized by chromosomal concentration, DNA degradation and cell shrinkage (24). Apoptosis acts as a quality control mechanism for the maintenance of tissue homeostasis by eliminating defective cells (25). Apoptosis not only exists in certain physiological processes, but is also involved in numerous pathological degenerative diseases, including neurodegeneration and IDD (26,27). Boos et al (28) reported that the initiation of IDD is associated with the changes of intervertebral disc cell behavior, including increased cell death. Gruber and Hanley (29) suggested that a large proportion of intervertebral disc cells underwent apoptosis in patients with degenerated intervertebral discs. Therefore, pathways of apoptosis may be involved in the progression of IDD in the ENP.
In addition, the present study found that the TF NFKB1 was upregulated in the pathway of apoptosis. In particular, NFKB1 was a hub gene in the PPI network. NFKB1 encodes the nuclear factor kappa B p105/p50 isoforms (30). NFKB proteins are a family of transcription factors that have critical roles in the regulation of various biological defense processes, including immune responses, cell-growth control and apoptosis (31). Inappropriate activation of NFKB has been associated with numerous inflammatory diseases (30). Wuertz et al (32) have reported that IDD is characterized not only by an imbalance between anabolic and catabolic processes, but also by inflammatory mechanisms. Of note, IL-1β and TNF-α are pro-inflammatory cytokines, which have been detected in the degenerated disc (33). In conclusion, the upregulation of NFKB1 may have been stimulated by IL-1β and TNF-α, which suggests that the TF NFKB1 may be an important biomarker for IDD.
In the present study, several downregulated DEGs, including COL1A1, were enriched in the pathway of ECM - receptor interaction; furthermore, COL1A1 was a hub gene in the PPI network. Over the previous decade, IDD research has focused on elucidating the mechanisms of ECM degradation, as it causes marked structural changes, including dehydration and fibrosis of the nucleus pulposus, disorganization of the annulus fibrosus and calcification of the cartilaginous end plates (28). These changes ultimately lead to structural failure in IDD. Of note, the protein collagen type I, alpha 1, encoded by COL1A1, is a major ECM component (34). Col1a1 has been found in most connective tissues and is abundant in bone (35). Several studies have demonstrated that the expression of genes encoding connective tissue components, particularly collagen, can be affected by numerous cytokines, including IL-1 and TNF-α (36,37). Mori et al (38) found that TNF-α decreased Col1a1 expression through suppressing promoter activity of Col1a1. Feng et al (39) suggested that Col1a1 can provide tensile strength in the annulus fibrosus of the intervertebral disc. Thus, the downregulated expression of COL1A1 in the present study may be regulated by IL-1β and TNF-α, resulting in the initiation of IDD.
In the PPI network, IL6 and HIF1A were hub genes with a high connectivity degree, while HIF1A was also a TF. IL6 together with IL-1α, IL-1β and TNF-α are among the most potent catabolic cytokines and pro-inflammatory mediators as mentioned above (40). Noponen-Hietala et al (41) reported that the features of pain, tissue destruction and inflammation in IDD were linked with the functions of IL6, and IL6 was also associated with IDD-associated radiculopathy. Therefore, IL6 is not the sole pro-inflammatory cytokine in the pathogenesis of IDD and appears to be important in the mediation of pain in IDD.
In the present study, the TF of HIF1A was also found in module 2 and was involved in GO terms associated with cell proliferation. HIF1 functions as a transcriptional regulator of the adaptive response to hypoxia and activates transcription of certain genes involved in energy metabolism, angiogenesis and apoptosis (42). Wu et al (43) found that Hif1a has an important role in the metabolism and synthesis of ECM and is a pivotal contributor to the survival of nucleus pulposus cells. Therefore, HIF1A may be involved in the development of IDD.
In conclusion, the present study provided a comprehensive bioinformatics analysis of the molecular mechanisms of IDD induced by IL-1β and TNF-α. Pathways of apoptosis and ECM - receptor interaction as well as their enriched DEGs may have important roles in the development and progression of IDD. In addition, the TFs NFKB1 and HIF1A and the DEGs COL1A1 and IL6 are hypothesized to interact with IL1-β and TNF-α to participate in the development of IDD. However, further bioinformatics and experimental studies with larger sample sizes are required to confirm the results of the present study.