This article is mentioned in:

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.

Table I GO functional enrichment analysis for the upregulated and downregulated differentially expressed genes.

Table I

GO functional enrichment analysis for the upregulated and downregulated differentially expressed genes.

GO Term	Biological process	Count	P-value
Enrichment score, 7.038
GO:0009611	Response to wounding	24	2.10×10−8
GO:0006954	Inflammatory response	17	5.50×10−8
GO:0006952	Defense response	21	6.67×10−7
Enrichment score, 5.439
GO:0010033	Response to organic substance	34	4.24×10−7
GO:0009719	Response to endogenous stimulus	24	4.27×10−6
GO:0009725	Response to hormone stimulus	21	2.66×10−5
Enrichment score, 6.414
GO:0007155	Cell adhesion	28	7.91×10−9
GO:0022610	Biological adhesion	28	7.91×10−9
GO:0016337	Cell-cell adhesion	12	9.22×10−4
Enrichment score, 6.331
GO:0030199	Collagen fibril organization	8	3.80×10−8
GO:0030198	Extracellular matrix organization	12	1.74×10−7
GO:0043062	Extracellular structure organization	12	1.54×10−5

[i] GO, gene ontology.

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.

Table II KEGG pathway enrichment analysis for the upregulated and downregulated DEGs.

Table II

KEGG pathway enrichment analysis for the upregulated and downregulated DEGs.

Term	Biological process/pathway	Count	P-value
Upregulated DEGs
rno04210	Apoptosis	9	4.15×10−5
rno04621	NOD-like receptor signaling pathway	7	3.43×10−4
rno04060	Cytokine-cytokine receptor interaction	11	8.06×10−4
rno04062	Chemokine signaling pathway	9	4.53×10−3
rno04630	Jak-STAT signaling pathway	8	5.29×10−3
rno05222	Small cell lung cancer	6	8.70×10−3
rno05200	Pathways in cancer	12	8.81×10−3
rno04620	Toll-like receptor signaling pathway	6	1.21×10−2
rno00230	Purine metabolism	7	3.06×10−2
Downregulated DEGs
rno04512	Extracellular matrix - receptor interaction	16	7.56×10−12
rno04510	Focal adhesion	20	1.13×10−9
rno05200	Pathways in cancer	14	4.34×10−3
rno04640	Hematopoietic cell lineage	6	1.30×10−2
rno04670	Leukocyte transendothelial migration	7	1.73×10−2
rno05410	Hypertrophic cardiomyopathy	6	1.75×10−2
rno04110	Cell cycle	7	2.59×10−2
rno05219	Bladder cancer	4	2.66×10−2
rno04115	p53 signaling pathway	5	3.08×10−2

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; DEG, differentially expressed gene.

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

Figure 1 Protein-protein interaction network of DEGs. Circular nodes indicate upregulated DEGs and square nodes signify downregulated DEGs. The color intensity is inversely proportional to the P-value. DEG, differentially expressed gene.

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.

Figure 2 Protein-protein interaction network of module 2. Circles and squares indicate DEGs and triangles signify transcription factors (upregulated). Circular nodes indicate upregulated DEGs and square nodes signify downregulated DEG. Arrows indicate transcriptional regulation, while lines signify protein-protein interactions. DEG, differentially expressed gene.

Table III GO and KEGG pathway enrichment analyses for the differentially expressed genes in module 2.

Table III

GO and KEGG pathway enrichment analyses for the differentially expressed genes in module 2.

Term	Biological process/pathway	Count	P-value
GO
GO:0042127	Regulation of cell proliferation	12	6.47×10−10
GO:0009611	Response to wounding	10	4.98×10−9
GO:0042060	Wound healing	7	1.89×10−7
GO:0001666	Response to hypoxia	7	3.28×10−7
GO:0070482	Response to oxygen levels	7	4.70×10−7
GO:0009725	Response to hormone stimulus	9	4.75×10−7
GO:0031099	Regeneration	6	8.12×10−7
GO:0009719	Response to endogenous stimulus	9	1.15×10−6
GO:0008284	Positive regulation of cell proliferation	8	1.20×10−6
GO:0001817	Regulation of cytokine production	6	3.26×10−6
KEGG
rno05200	Pathways in cancer	8	7.66×10−6
rno05212	Pancreatic cancer	4	7.37×10−5
rno04510	Focal adhesion	5	1.45×10−3
rno04150	mTOR signaling pathway	3	8.56×10−3
rno04621	NOD-like receptor signaling pathway	3	1.16×10−2
rno04115	p53 signaling pathway	3	1.31×10−2
rno05211	Renal cell carcinoma	3	1.42×10−2
rno04060	Cytokine-cytokine receptor interaction	4	1.43×10−2
rno05220	Chronic myeloid leukemia	3	1.67×10−2
rno05222	Small cell lung cancer	3	2.02×10−2
rno05215	Prostate cancer	3	2.35×10−2

[i] KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, gene ontology.

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.

Acknowledgments

The present study was supported by the National Natural Youth Science Foundation of China (grant no. 31300778).

References