Identifying the potential role of IL-1β in the molecular mechanisms of disc degeneration using gene expression profiling and bioinformatics analysis

Abstract

Purpose

We performed a bioinformatics analysis to identify the key genes that were differentially expressed between degenerative intervertebral disc (IVD) cells with and without exposure to interleukin-1β and explore the related signaling pathways and interaction networks.

Methods

The microarray data were downloaded from the Gene Expression Omnibus (27,494). Then, analyses of the gene ontology, signaling pathways, and interaction networks for the differentially expressed genes (DEGs) were conducted using tools including the Database for Annotation, Visualization, and Integrated Discovery, Metascape, Gene Set Enrichment Analysis, Search Tool for the Retrieval of Interacting Genes, Cytoscape, Venn method, and packages of the R computing language.

Results

A total of 260 DEGs were identified, including 161 upregulated and 99 downregulated genes. Gene Ontology annotation analysis showed that these DEGs were mainly associated with the extracellular region, chemotaxis, taxis, cytokine activity, and cytokine receptor binding. A Kyoto Encyclopedia of Genes and Genomes signaling pathway analysis showed that these DEGs were mainly involved in the of cytokine-cytokine receptor interaction, rheumatoid arthritis, tumor necrosis factor signaling pathway, Salmonella infection, and chemokine signaling pathway. The interaction network analysis indicated that 10 hub genes, including CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, PF4, and GPER1 may play key roles in IVD degeneration.

Conclusions

Bioinformatic analysis showed that CXCL8 and other nine key genes may play a role in the development of disc degeneration induced by inflammatory reactions and can be used to identify potential target genes for therapeutic applications in IVD degeneration.

Keywords

intervertebral disc degeneration interleukin-1β gene expression profiling computational biology

Introduction

Intervertebral disc degeneration (IVDD) is considered to be the leading cause of low back pain, severely threatening the quality of life of middle-aged and elderly people and increasing the economic burden on families and the society.^1,2 Although modern evidence-based medicine has identified IVDD as a consequence of a variety of genetic, traumatic, lifestyle, aging, and nutritional factors, the pathological mechanisms involved in IVDD progression remain unknown.^3–7 Currently, treatment modalities include noninvasive treatments, such as drug therapy, multiple physical modalities, and multidisciplinary biopsychosocial rehabilitation; interventional modalities, such as intradiscal therapies and epidural injections; regenerative modalities with disc injections of various solutions; and surgical approaches, such as fusion and artificial disc replacement. Despite the progress made in therapies aimed at relieving the pain symptoms, they offer only temporary benefits and are related to complications.⁸ Therefore, it is necessary to elucidate the underlying pathogenesis of IVDD and lower back pain to develop novel effective treatment methods for these conditions.

The role of bioinformatics in the exploration of disease mechanisms has become essential given the multifactorial nature of many common health conditions. Examples include the identification of hub genes associated with the etiology of intervertebral disc (IVD).^9–12 Inflammatory cytokines play a significant role in disc degeneration.¹³ As disc degeneration progresses, the expression of pro-inflammatory cytokines, such as interleukin (IL)-1β and tumor necrosis factor alpha (TNF-α) in the nucleus pulposus (NP) tissue increases significantly, which could accelerate disc cell degeneration by promoting the release of matrix metalloproteinases (MMPs) and inflammatory cytokines.^13–16 However, Le et al. believed that IL-1 may play a more significant role than TNF-α.¹⁶ Furthermore, IL-1β has been shown to elevate both neurotrophin gene expression levels and the production of nerve growth factor, which is associated with low back pain.¹⁷ Thereafter, inflammatory processes exacerbated by IL-1β are believed to be key mediators of disc degeneration and low back pain, and IL-1β may be a better target for therapeutic intervention of IVDD.

In this study, microarray data profiling of human disc tissue samples from patients with herniated discs and degenerative disc disease, including four samples stimulated with 10² pM IL-1β and four controls performed by Gruber et al. (GSE27494) were downloaded and reanalyzed. Differentially expressed genes (DEGs) were identified and co-expression network analysis was performed, followed by functional analysis of the genes in the network. Differentially expressed genes were subsequently subjected to functional and pathway annotations, construction of protein-protein interaction (PPI) networks, and module analysis to investigate the critical DEGs in the progression of IVDD. The present study aimed to identify the molecular markers indicating the induction of disc degeneration by inflammatory cytokines, which may guide further investigations and contribute to the development of novel therapeutic strategies for the treatment of IVDD.

Materials and methods

Collection of microarray data

The gene profile GSE27494, which was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), was produced using the [U133_X3P] Affymetrix Human X3P Array (Platform GPL1352). The gene profile consisted of four degenerated disc tissue samples and four degenerative disc tissues subjected to additional IL-1β exposure. Disc tissue samples were obtained from surgical disc procedures performed on patients with herniated discs and degenerative disc disease. Cultured annulus fibrosus (AF) cells were grown in a three-dimensional collagen construct with or without 10² pM IL-1 for 14 days. Following homogenization in TRIzol reagent, total RNA was isolated and analyzed using microarrays.

Identification of DEGs

The limma R package has emerged as one of the most widely used statistical tests to identify DEGs.¹⁸ The package was obtained from an open website (http://www.bioconductor.org/packages/release/bioc/html/limma.html). We screened DEGs between IL-1β-exposed and control disc cells by utilizing the limma package with an adjusted p-value < 0.05 and a log (fold change) > 2 or log (fold change) < −2 as the cut-off criteria. Furthermore, volcano plots and heat maps were generated using the R software.

Functional annotation of DEGs using the database for annotation, visualization, and integrated discovery, and Metascape analysis

We used Database for Annotation, Visualization, and Integrated discovery (DAVID) (https://david.ncifcrf.gov/home.jsp) (version 6.8)¹⁹ for Gene Ontology (GO) functional analysis²⁰ and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs.²¹ The enriched GO terms and pathways with a p-value < 0.05 were identified. Moreover, GO and KEGG pathway analyses were performed using Metascape (http://metascape.org/gp/index.html#/main/step1).²² Terms with a p-value < 0.01, a minimum count of 3, and an enrichment factor >1.5 (enrichment factor is the ratio between the observed count and the count expected by chance) were collected and grouped into clusters based on their membership similarities.

Enrichment analysis of GO and KEGG using gene set enrichment analysis

Gene set enrichment analysis primarily analyzes microarray data using genomic and genetic sequencing to detect significant biological differences in microarray datasets.²³ In the present study, we performed GSEA analysis on the gene sequences of IL-1β-exposed and control disc cells using c5.all.v7.1.symbols.gmt, and c2.cp.kegg.v7.1.symbols.gmt, respectively, as reference gene sets (software.broadinstitute.org/gsea/).²⁴ Gene set permutations were analyzed 1000 times. In addition, the normalized enrichment score, normalized p-value, and false discovery rate q values were used to filter the correlative GO terms and pathways.

Construction of the PPI network, significant module, and hub gene network

First, we used Metascape (http://metascape.org/gp/index.html#/main/step1)²² to construct the PPI network and screen the significant module. Second, to assess the functional associations among the target genes of upregulated and downregulated DEGs, we uploaded the target gene data to the Search Tool for the Retrieval of Interacting Genes (STRING, http://string.embl.de/) to construct the PPI network. Interactions with a combined score >0.7 were considered significant. Highly interconnected (hub) genes in the PPI network were analyzed using Cytoscape (version 3.6.0). The Molecular Complex Detection tool (MCODE) (version 1.5.1),²⁵ a plug-in of Cytoscape, was used to screen and identify the most significant module in the PPI network. Genes with the top 10 highest degrees in the PPI network were identified as the hub genes.

Identification of significant genes

We generated a Venn diagram to identify the significant common genes among “Metascape_MCODE,” “Cytoscape_MCODE,” and “Cytoscape_cytoHubba” using Fun Rich software (http://www.funrich.org).²⁵ Summaries of the functions of the significant genes were obtained via DAVID (https://david.ncifcrf.gov/home.jsp). The R language was used to perform the clustering analysis of significant genes based on the gene expression levels.

Results

Screening of DEGs

After the analysis of the datasets (GSE27494) using the limma package, the DEGs between IL-1β-exposed and control disc cells were presented in a volcano plot (Figure 1). After setting the threshold value, a total of 260 DEGs were reserved, including 161 upregulated DEGs (Table 1) and 99 downregulated DEGs (Table 2). The cluster heatmap of the DEGs is shown in Figure 2.

Figure 1.

Volcano plot presenting the DEGs between IL-1β-exposed and control disc cells. The downregulated DEGs are marked in green, and the upregulated DEGs are marked in red. DEG: differentially expressed genes.

Table 1.

The genes that were upregulated by a fold change > 2 and with a p < 0.05 (top 20).

Gene symbol	Ensembl gene ID	Log FC	p-value	adj.P.Val
CXCL5	ENSG00000163735	9.83	9.95 × 10–06	0.017479
CXCL8	ENSG00000169429	7.88	4.72 × 10–06	0.012986
IL1B	ENSG00000125538	7.29	3.10 × 10–05	0.027488
MMP3	ENSG00000149968	7.18	1.68 × 10–04	0.057873
SLC7A2	ENSG00000003989	6.90	4.16 × 10–06	0.012986
C15orf48	ENSG00000166920	6.27	1.41 × 10–05	0.021941
CXCL3	ENSG00000163734	6.24	1.02 × 10–05	0.017479
CCL3L3	ENSG00000276085	6.20	1.04 × 10–03	0.120405
CXCL2	ENSG00000081041	6.18	1.62 × 10–05	0.022413
IL11	ENSG00000095752	6.10	5.88 × 10–05	0.036427
IL6	ENSG00000136244	5.87	2.27 × 10–03	0.169458
HSD11B1	ENSG00000117594	5.49	4.31 × 10–05	0.03215
LIF	ENSG00000128342	5.35	8.59 × 10–06	0.016999
RSPO3	ENSG00000146374	5.35	2.46 × 10–05	0.027393
MMP8	ENSG00000118113	5.29	4.37 × 10–06	0.012986
CLDN1	ENSG00000163347	5.17	9.33 × 10–06	0.017479
ERC2	ENSG00000187672	5.11	3.09 × 10–05	0.027488
EHF	ENSG00000135373	5.02	2.05 × 10–04	0.063472
BMP2	ENSG00000125845	4.89	6.30 × 10–04	0.098134
RAB27B	ENSG00000041353	4.86	7.02 × 10–06	0.01485

FC: fold change; MMP: matrix metalloproteinases.

Table 2.

The genes that were downregulated by a fold change < −2 and with a p < 0.05 (top 20).

Gene symbol	Ensembl gene ID	Log FC	p-value	adj.P.Val
FGL2	ENSG00000127951	−5.42	9.61 × 10–05	0.047186
LSP1	ENSG00000130592	−5.03	6.50 × 10–04	0.100144
ACTC1	ENSG00000159251	−4.94	2.86 × 10–02	0.48192
FAM107A	ENSG00000168309	−4.75	9.74 × 10–03	0.307323
MAP1LC3C	ENSG00000197769	−4.67	4.10 × 10–04	0.080189
CRIP1	ENSG00000213145	−4.59	7.58 × 10–04	0.107444
EXOSC7	ENSG00000075914	−4.49	2.16 × 10–05	0.025494
COL15A1	ENSG00000204291	−4.34	7.55 × 10–03	0.275959
COMP	ENSG00000105664	−4.29	4.53 × 10–02	0.573519
ASPN	ENSG00000106819	−4.10	1.03 × 10–02	0.316056
ADIRF	ENSG00000148671	−4.03	1.27 × 10–03	0.132417
OGN	ENSG00000106809	−4.00	1.51 × 10–03	0.14488
SMOC2	ENSG00000112562	−3.96	2.70 × 10–04	0.070341
OLFML2B	ENSG00000162745	−3.79	1.75 × 10–04	0.058275
PDE5A	ENSG00000138735	−3.76	3.24 × 10–04	0.072076
SLC14A1	ENSG00000141469	−3.63	4.25 × 10–02	0.557619
SNED1	ENSG00000162804	−3.56	1.72 × 10–03	0.152379
KCNN4	ENSG00000104783	−3.35	2.59 × 10–03	0.176067
TOX	ENSG00000198846	−3.34	5.00 × 10–03	0.23011
ANK3	ENSG00000151150	−3.30	2.19 × 10–03	0.167475

FC, fold change.

Figure 2.

Heatmap of DEGs between IL-1β-exposed and control disc cells. Each column represents a DEG and each row represents a single sample. Red color indicates an upregulation and the green color indicates a down-regulation. DEG: differentially expressed genes.

The enrichment results of GO and KEGG via DAVID and Metascape

The GO functional annotation analysis of the DEGs showed that (1) the biological processes were mainly involved in the regulation of chemotaxis, taxis, and the response to molecules of bacterial origin; (2) the molecular functions of the altered genes were mainly related to cytokine activity, cytokine receptor binding, and chemokine receptor binding; and (3) the cell components involved were mainly in the extracellular region and extracellular matrix (Figure 3(a)–(c)). In the KEGG signaling pathway analysis, the 260 DEGs were found to be mainly involved in cytokine-cytokine receptor interaction, rheumatoid arthritis, TNF signaling pathway, Salmonella infection, chemokine signaling pathway, and other signal transduction pathways (Figure 3(d)).

Figure 3.

Bubble diagram of GO analysis and KEGG pathway analysis using DAVID. Detailed information relating to changes in the (a) biological process, (b) molecular function, (c) cellular components, and (d) KEGG analysis for hub genes. KEGG: Kyoto Encyclopedia of Genes and Genomes; DAVID: Database for Annotation, Visualization and Integrated discovery; GO: Gene Ontology.

Furthermore, the functional enrichment analysis performed using Metascape showed that the DEGs were significantly enriched in NABA MATRISOME ASSOCIATED, extracellular matrix organization, NABA CORE MATRISOME, positive regulation of cell motility, and regulation of cell adhesion (Figure 4(a)–(c)).

Figure 4.

Enrichment analysis of differentially expressed genes using Metascape. (a) Heatmap of enriched terms across input differentially expressed gene lists, colored by p-values, obtained using Metascape. (b) Network of enriched terms colored by cluster identity, in which the nodes that share the same cluster identity are typically close to each other. (c) Network of enriched terms colored by p-value, in which the terms containing more genes tend to have a more significant p-value.

Gene Ontology and KEGG pathway enrichment analysis of DEGs using GSEA

The results of GO enrichment analysis using GSEA indicated that 3186/5154 genes were upregulated in IL-1β-exposed disc cells. The most significantly enriched genes in the upregulated and downregulated gene sets according to the normalized enrichment score are listed in Table 3. Figure 5 shows the six most significant plots in the upregulated and downregulated gene sets. Kyoto Encyclopedia of Genes and Genomes enrichment analysis using GSEA indicated that 107/177 gene sets were upregulated in IL-1β-exposed disc cells compared to control disc cells, while 70/177 gene sets were downregulated. The most significantly enriched genes in the upregulated and downregulated gene sets according to the normalized enrichment score are listed in Table 4. Figure 6 shows the six most significant plots in the upregulated and downregulated gene sets.

Table 3.

Functional enrichment analysis of DEGs in IL-1β-exposed disc cells using GSEA.

Gene set name	Size	ES	NES	p-value
Upregulated
GO_NEUTROPHIL_MIGRATION	103	0.835	2.461	0.000
GO_ANTIMICROBIAL_HUMORAL_IMMUNE_RESPONSE_MEDIATED_BY_ANTIMICROBIAL_PEPTIDE	58	0.895	2.387	0.000
GO_GRANULOCYTE_MIGRATION	125	0.789	2.387	0.000
GO_CYTOKINE_ACTIVITY	194	0.750	2.366	0.000
GO_CHEMOKINE_RECEPTOR_BINDING	50	0.891	2.354	0.000
GO_CHEMOKINE_ACTIVITY	41	0.908	2.335	0.000
Downregulated
GO_REGULATION_OF_INOSITOL_PHOSPHATE_BIOSYNTHETIC_PROCESS	16	0.842	1.810	0.004
GO_IMMUNOGLOBULIN_RECEPTOR_BINDING	15	0.852	1.795	0.000
GO_WATER_TRANSPORT	16	0.821	1.753	0.000
GO_CEREBRAL_CORTEX_NEURON_DIFFERENTIATION	24	0.763	1.735	0.002
GO_EXTRACELLULAR_MATRIX_STRUCTURAL_CONSTITUENT	157	0.558	1.731	0.000
GO_REGULATION_OF_MEMBRANE_DEPOLARIZATION	43	0.676	1.728	0.000

ES: Enrichment Score; NES: Normalized Enrichment Score; DEG: differentially expressed genes; GSEA: Gene Set Enrichment Analysis.

Figure 5.

Six significant enrichment plots of functional enrichment analysis of DEGs between IL-1β-exposed and control samples in GSE27494 using GSEA. (a) Enrichment plot: GO_REGULATION_OF_INOSITOL_PHOSPHATE_BIOSYNTHETIC_PROCESS. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (b) Enrichment plot: GO_IMMUNOGLOBULIN_RECEPTOR_BINDING. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (c) Enrichment plot: GO_WATER_TRANSPORT. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (d) Enrichment plot: GO_NEUTROPHIL_MIGRATION. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (e) Enrichment plot: GO_ANTIMICROBIAL_HUMORAL_IMMUNE_RESPONSE_MEDIATED_BY_ANTIMICROBIAL_PEPTIDE. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (f) Enrichment plot: GO_GRANULOCYTE_MIGRATION. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. DEG: differentially expressed genes; GSEA: Gene Set Enrichment Analysis.

Table 4.

Pathway enrichment analysis of DEGs using GSEA.

Gene set name	Size	ES	NES	p-value
Upregulated
KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION	243	0.742	2.420	0.000
KEGG_CHEMOKINE_SIGNALING_PATHWAY	171	0.665	2.075	0.000
KEGG_GRAFT_VERSUS_HOST_DISEASE	35	0.823	2.035	0.000
KEGG_NOD_LIKE_RECEPTOR_SIGNALING_PATHWAY	62	0.740	2.002	0.000
KEGG_TOLL_LIKE_RECEPTOR_SIGNALING_PATHWAY	96	0.654	1.931	0.000
KEGG_HEMATOPOIETIC_CELL_LINEAGE	75	0.678	1.900	0.000
Downregulated
KEGG_BASE_EXCISION_REPAIR	33	0.617	1.493	0.033
KEGG_THYROID_CANCER	29	0.607	1.475	0.045
KEGG_DILATED_CARDIOMYOPATHY	85	0.498	1.422	0.021
KEGG_DNA_REPLICATION	36	0.536	1.351	0.069
KEGG_PEROXISOME	75	0.438	1.243	0.112
KEGG_TGF_BETA_SIGNALING_PATHWAY	84	0.423	1.219	0.113

ES: Enrichment Score; NES: Normalized Enrichment Score; DEG: differentially expressed genes; GSEA: Gene Set Enrichment Analysis.

Figure 6.

Six significant enrichment plots of pathway enrichment analysis of DEGs between IL-1β-exposed and control samples in GSE27494 using GSEA. (a) Enrichment plot: KEGG_THYROID_CANCER. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (b) Enrichment plot: KEGG_DILATED_CARDIOMYOPATHY. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (c) Enrichment plot: KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (d) Enrichment plot: KEGG_BASE_EXCISION_REPAIR. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (e) Enrichment plot: KEGG_CHEMOKINE_SIGNALING_PATHWAY. Profile of the Running ES Score & Positions of Gene Set Members on the Rank Ordered List. (f) Enrichment plot: KEGG_GRAFT_VERSUS_HOST_DISEASE. Profile of the Running ES Score & Positions of GeneSet Members on the Rank Ordered List. DEG: differentially expressed genes; GSEA: Gene Set Enrichment Analysis.

Protein-protein interaction construction and module analysis

Using Metascape analysis, a PPI network of DEGs was constructed (Figure 7(a)). Six MCODE modules were identified in the PPI network. MCODE 1 consisted of 12 genes, including CXCL2, CCL20, SAA1, G protein-coupled estrogen receptor 1 (GPER1), CXCL8, CXCL1, C3, ACKR3, PF4, CXCL5, CXCL3, and CXCL6 (Figure 7(b)). Molecular Complex Detection tool 2 consisted of five genes, including ACTC1, GCH1, TUBB3, ACTG2, and SOD2 (Figure 7(c)). Molecular Complex Detection tool 3 consisted of three genes, TF, TFRC, and AMPH (Figure 7(d)). Molecular Complex Detection tool 4 consisted of three genes: Lyn, IL6, and IL1B (Figure 7(e)). Molecular Complex Detection tool 5 consisted of three genes: MMP3, CCL7, and CCL2 (Figure 7(f)). Lastly, MCODE 6 consisted of three genes: CCKAR, GNA14, and P2RY6 (Figure 7(g)).

Figure 7.

PPI network and significant module constructed using Metascape. (a) Through Metascape analysis, a protein-protein interaction network of DEGs was constructed. (b) MCODE 1 consists of 12 genes: CXCL2, CCL20, SAA1, GPER1, CXCL8, CXCL1, C3, ACKR3, PF4, CXCL5, CXCL3, and CXCL6. (c) MCODE 2 consists of 5 genes: ACTC1, GCH1, TUBB3, ACTG2, and SOD2. (d) MCODE 3 consists of 3 genes: TF, TFRC, and AMPH. (e) MCODE 4 consists of 3 genes: LYN, IL6, and IL1B. (f) MCODE 5 consists of 3 genes: MMP3, CCL7, and CCL2. (g) MCODE 6 consists of 3 genes: CCKAR, GNA14, and P2RY6. PPI: protein-protein interaction; DEG: differentially expressed genes.

Furthermore, the PPI network of the DEGs was also generated using the STRING online tool, and 327 edges and 251 nodes were observed (PPI enrichment p-value < 0.05) (Figure 8(a)). Six MCODE modules were also identified from the PPI network using the MCODE (version 1.5.1), one plug-in of Cytoscape (Figure 8(b)–(g)). Using Cytoscape, a total of 10 genes were identified as hub genes: CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, PF4, and GPER1 (Figure 8(h)).

Figure 8.

PPI network, significant module, and hub gene network constructed using STRING and Cytoscape. (a) There were 327 edges and 251 nodes in the PPI network. All nodes with a combined score >0.7 are shown. Disconnected nodes in the network were excluded. (b) The first module consists of 95 edges and 17 nodes. (c) The second module consists of 15 edges and 6 nodes. (d) The third module consists of 10 edges and 5 nodes. (e) The fourth module consists of 3 edges and 3 nodes. (f) The fifth module consists of 3 edges and 3 nodes. (g) The sixth module consists of 3 edges and 3 nodes. (h) A total of 10 genes were identified as hub genes using Cytoscape: CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, PF4, and GPER1. PPI: protein-protein interaction; GPER1: G protein-coupled estrogen receptor 1; STRING: Search Tool for the Retrieval of Interacting Genes.

Identification and analysis of significant genes

The VENN diagram showed that there were 10 significant common genes among “Metascape_MCODE,” “Cytoscape_MCODE,” and “Cytoscape_cytoHubba,” including: CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, PF4, and GPER1 (Figure 9(a)). Hierarchical clustering showed that the significant genes could be used to significantly distinguish IL-1β-exposed disc cells from control cells. Compared with the control cells, the expression of GPER1 was downregulated, and the expression of CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, and PF4 were upregulated in IL-1β-exposed discs (Figure 9(b)).

Figure 9.

Identification and analysis of significant genes. (a) The VENN diagram shows that there were 10 significant common genes among “Metascape_MCODE,” “Cytoscape_MCODE,” and “Cytoscape_cytoHubba,”: CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, PF4, and G protein-coupled estrogen receptor 1. (b) Hierarchical clustering showed that the significant genes could basically differentiate the IL-1β-exposed disc cells from the control disc cells.

Discussion

Biological information analysis has been widely applied to help scientists identify new therapeutic targets by exploring gene changes in the process of disease development.^26,27 In the present study, we focused on inflammation as the primary etiological factor impacting disc degeneration, and IL-1β was selected as a representative inflammatory cytokine.²⁸ Therefore, we obtained gene expression profiling data from GEO microarray datasets containing four samples of degenerative IVD cells stimulated with 10² pM IL-1β and four controls. The key genes that were affected during disc degeneration caused by inflammatory cytokine exposure were analyzed using bioinformatics analysis, and 260 DEGs were obtained (161 upregulated and 99 downregulated). Subsequently, GO annotation and KEGG pathway enrichment analyses of the DEGs revealed their participation in several cellular processes (e.g., the regulation of chemotaxis and taxis) and signaling cascades related to IVDD development (e.g., the interactions of cytokine-cytokine receptor interaction and rheumatoid arthritis). A PPI network based on these genes was constructed to obtain the top 10 hub genes: CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, PF4, and GPER1. Our results suggest that increased levels of IL-1β affected the expression of the hub genes, and subsequent analysis confirmed their involvement in important IVDD-related pathways. Therefore, we believe that these hub genes might be potential biological targets for IVDD diagnosis and for the development of therapeutic drugs.

Chemokines produced by mammalian cells during inflammation are a class of chemotactic and inducible small molecule peptides that are ubiquitous and play an important role in acute and chronic inflammation.²⁹ Based on key cysteine residues involved in disulfide bonds, chemokines are classified as C, CC, CXC, and CX3C.³⁰ Among the top 10 hub genes observed in the present study, CXCL8, CXCL1, CXCL2, CXCL5, CXCL3, CXCL6, and PF4 belong to the CXC family of chemokines, CCL20 belongs to the CC family, and C3 belongs to the C family. Of the nine chemokines, CXCL8, CXCL6, and CCL20 have been confirmed to participate in IVDD in previous studies. Zhang et al.³¹ found that AF samples from patients with low back pain had an elevated IL-8 expression compared to controls with scoliosis. CXCL6 has also been detected in a conditioned medium of induced-degenerative discs in an organ culture and may contribute to the chemotactic response of induced-degenerative discs.³² Grad et al.³³ observed significantly elevated levels of systemic blood plasma concentrations of chemokine CXCL6 in individuals with moderate/severe disc degeneration, in comparison with a control group with no signs of disc degeneration, according to the MRI results. Meanwhile, Zhang et al.³⁴ reported that NP cells from the extruded and herniated patient groups could produce abundant amounts of CCL20. Subsequently, Zhang et al.³⁵ suggested a potential mechanism for the recruitment of IL-17-producing cells to degenerated IVDs via a CCL20/CCR6 system in vivo, while IL-17 is involved in the autoimmune process of IVDD in rat models. In the present study, the expression of CXCL8, CXCL6, and CCL20 in IL-1β-exposed samples were increased by 7.88-fold, 3.2-fold, and 4.18-fold, respectively, compared to control samples, respectively. Therefore, we speculated that increased levels of IL-1β affected the expression of CXCL8, CXCL6, and CCL20 in IVD cells, suggesting that these chemokines may be potential new targets for treating IVDD. However, we did not retrieve any literature identifying whether IVDD is correlated with CXCL5, CXCL1, CXCL2, CXCL3, C3, and PF4. Therefore, further research investigating the six genes is needed to determine their detailed mechanisms in IVDD.

Previous studies have shown that pathological inflammation of IVD and peridiscal space is characterized by increased levels of pro-inflammatory cytokines and chemokines.^36–38 Cytokines and chemokines have three modes of action: (1) recruiting inflammatory cells and activation of phagocytosis, (2) stimulating the production of other inflammatory mediators and MMPs, and (3) enhancing matrix degradation.^39–44 In terms of GO enrichment analysis, we found that most DEGs were mainly involved in the regulation of chemotaxis and taxis, cytokine activity, cytokine receptor binding, and chemokine receptor binding. Regarding the KEGG pathway, we found that most of the DEGs were primarily involved in cytokine-cytokine receptor interaction, TNF signaling pathway, and chemokine signaling pathway. In addition, among the top 10 hub genes, 9 chemokines were the most significantly upregulated genes. Therefore, in addition to the direct effect on IVDD, increased levels of IL-1β may activate a cascade of detrimental self-promoting events that exacerbate IVDD by stimulating the production of pro-inflammatory cytokines and chemokines. Furthermore, experimental studies have reported the role of pro-inflammatory cytokines in degenerated and painful IVD,^45,46 which indicates that IL-1β is significantly associated with low back pain.

The underlying molecular mechanisms of IVDD have been previously investigated, and extracellular matrix degradation and inflammation have been confirmed to play a critical role in accelerating IVDD progression. A normal adult IVD is composed of vascular tissue, which contains a large amount of extracellular matrix, and the disc cells are responsible for maintaining the integrity of the extracellular matrix.⁴⁷ Changes in the metabolic balance of IVD cells affect the quality and quantity of the extracellular matrix and its functional properties; therefore, these changes can be related to disc degeneration.⁴⁸ Our study showed that the DEGs were mostly enriched in the extracellular region and extracellular matrix, which indicates that IL-1β may be involved in IVDD by affecting the quality and quantity of the extracellular matrix. In addition, Le et al.¹⁶ reported that IL-1β induced the expression of MMPs, and the upregulation of MMPs expression is implicated in disc extracellular matrix destruction, leading to the development of IDD.⁴⁹ Similarly, our bioinformatics analysis showed that MMP-3 significantly increased in degenerative IVD cells stimulated with 10² pM IL-1β compared to controls. Therefore, we inferred that IL-1β may affect the extracellular matrix and participate in IVDD by stimulating the expression of MMPs, although MMP-3 is not among the top 10 hub genes.

Furthermore, GPER1, which is an estrogen receptor (ER), was the only downregulated hub gene (downregulated by 2.42-fold). There are four major naturally occurring estrogens identified in women: oestrone (E1), estradiol (E2), oestriol (E3), and estetrol (E4). In addition, ERs consist of two classic nuclear receptors, ER-α and ER-β, and membrane-bound GPR30.⁵⁰ Over the past decade, many researchers have studied the relationship between estrogen and IVDD and found that estrogen can effectively alleviate IVDD development by inhibiting the apoptosis of IVD cells. Eestrogen can decrease IVD cell apoptosis in multiple ways, including the inhibition of inflammatory cytokines IL-1β and TNF-α, reducing catabolism due to MMP inhibition, upregulating integrin α2β1 and IVD anabolism, activating the PI3K/Akt pathway, decreasing oxidative damage, and promoting autophagy.⁵⁰ Song et al. found that both cytoplasmic and nuclear staining of ER-α and ER-β immunoreactivity were observed in NP cells, and ER-α and ER-β expression significantly decreased, along with the aggravation of IVDD in both males and females.⁵¹ Wei et al.⁵² also showed that GPR30, which has a high affinity for estrogen, is expressed in the human disc NP and can mediate E2-enhanced cell proliferation and influence disc cell survival. Therefore, we speculate that IL-1β decreases the expression of GPER1 in human disc AF, which accelerates disc degeneration.

Despite the rigorous bioinformatic analysis performed in the present study, there were still some limitations. First, the size of the dataset was small; thus, further analyses with an increased sample size are necessary to confirm the results. Second, we did not conduct further experimental verification of the screened genes. Although our study had these limitations, we believe that the results of the present study may provide potential directions for future research.

Conclusions

In this study, bioinformatic analysis of the gene expression profiles of degenerative IVD cells stimulated with IL-1β showed that CXCL8 and other nine key genes may play a role in the development of disc degeneration induced by inflammatory reactions. This suggests that bioinformatics methods can be used to identify potential therapeutic target genes and provide new insights into IVDD.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the key medical disciplines/schools of Shijingshan district, Beijing.

ORCID iD

Lei Zang

References

Vergroesen

Kingma

Emanuel

, et al. Mechanics and biology in intervertebral disc degeneration: a vicious circle. Osteoarthr Cartil 2015; 23(7): 1057–1070.

Liao

Song

, et al. Angiopoietin‐like protein 8 expression and association with extracellular matrix metabolism and inflammation during intervertebral disc degeneration. J Cell Mol Med 2019; 23(8): 5737–5750.

Teraguchi

Yoshimura

Hashizume

, et al. Progression, incidence, and risk factors for intervertebral disc degeneration in a longitudinal population-based cohort: the Wakayama Spine Study. Osteoarthr Cartil 2017; 25(7): 1122–1131.

Dario

Ferreira

Refshauge

, et al. The relationship between obesity, low back pain, and lumbar disc degeneration when genetics and the environment are considered: a systematic review of twin studies. Spine J 2015; 15(5): 1106–1117.

Sivan

Wachtel

Roughley

. Structure, function, aging and turnover of aggrecan in the intervertebral disc. Biochim Biophys Acta 2014; 1840(10): 3181–3189.

Hangai

Kaneoka

Kuno

, et al. Factors associated with lumbar intervertebral disc degeneration in the elderly. Spine J 2008; 8(5): 732–740.

Wang

YXJ

. Postmenopausal Chinese women show accelerated lumbar disc degeneration compared with Chinese men. J Orthop Translat 2015; 3(4): 205–211.

Zhao

Manchikanti

Kaye

, et al. Treatment of discogenic low back pain: current treatment strategies and future options-a literature review. Curr Pain Headache Rep 2019; 23(11): 86.

Sampara

Banala

Vemuri

, et al. Understanding the molecular biology of intervertebral disc degeneration and potential gene therapy strategies for regeneration: a review. Gene Ther 2018; 25(2): 67–82.

10.

Zhang

Chen

, et al. Melatonin modulates IL-1β-induced extracellular matrix remodeling in human nucleus pulposus cells and attenuates rat intervertebral disc degeneration and inflammation. Aging 2019; 11(22): 10499–10512.

11.

Wang

Chen

Ying

, et al. Ligustilide alleviated IL-1β induced apoptosis and extracellular matrix degradation of nucleus pulposus cells and attenuates intervertebral disc degeneration in vivo. Int Immunopharmacol 2019; 69: 398–407.

12.

Gorth

Shapiro

Risbud

. A new understanding of the role of IL‐1 in age‐related intervertebral disc degeneration in a murine model. J Bone Miner Res 2019; 34(8): 1531–1542.

13.

Risbud

Shapiro

. Role of cytokines in intervertebral disc degeneration: pain and disc content. Nat Rev Rheumatol 2014; 10(1): 44–56.

14.

Zhang

, et al. Simvastatin inhibits IL-1β-induced apoptosis and extracellular matrix degradation by suppressing the NF-kB and MAPK pathways in Nucleus Pulposus Cells. Inflammation 2017; 40(3): 725–734.

15.

Huang

Hsu

Chen

, et al. The roles of IL-19 and IL-20 in the inflammation of degenerative lumbar spondylolisthesis. J Inflamm 2018; 15: 19.

16.

Le Maitre

Hoyland

Freemont

. Catabolic cytokine expression in degenerate and herniated human intervertebral discs: IL-1β and TNFα expression profile. Arthritis Res Ther 2007; 9(4): R77.

17.

Gruber

Hoelscher

Bethea

, et al. Interleukin 1-beta upregulates brain-derived neurotrophic factor, neurotrophin 3 and neuropilin 2 gene expression and NGF production in annulus cells. Biotech Histochem 2012; 87(8): 506–511.

18.

Ritchie

Phipson

, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015; 43(7): e47.

19.

Huang

Sherman

Tan

, et al. The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol 2007; 8(9): R183.

20.

Mootha

Lindgren

Eriksson

, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genetics 2003; 34(3): 267–273.

21.

Szklarczyk

Franceschini

Wyder

, et al. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res 2015; 43(Database issue): D447–D452.

22.

Zhou

Pache

, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Communicat 2019; 10(1): 1523.

23.

Lai

Zhang

Nayak

, et al. Concordant integrative gene set enrichment analysis of multiple large-scale two-sample expression data sets. BMC Genomics 2014; 15 Suppl 1(Suppl 1): S6.

24.

Subramanian

Tamayo

Mootha

, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA 2005; 102(43): 15545–15550.

25.

Zhao

Quan

, et al. Identification of significant gene biomarkers of low back pain caused by changes in the osmotic pressure of nucleus pulposus cells. Scient Rep 2020; 10(1): 3708.

26.

Emmert-Streib

Zhang

Hamilton

. Computational cancer biology: education is a natural key to many locks. BMC Cancer 2015; 15: 7.

27.

Ferté

Trister

Huang

, et al. Impact of bioinformatic procedures in the development and translation of high-throughput molecular classifiers in oncology. Clin Cancer Res 2013; 19(16): 4315–4325.

28.

Phillips

Cullen

Chiverton

, et al. Potential roles of cytokines and chemokines in human intervertebral disc degeneration: interleukin-1 is a master regulator of catabolic processes. Osteoarthrit Cartil 2015; 23(7): 1165–1177.

29.

White

Iqbal

Greaves

. CC chemokine receptors and chronic inflammation--therapeutic opportunities and pharmacological challenges. Pharmacol Rev 2013; 65(1): 47–89.

30.

Beider

Abraham

Peled

. Chemokines and chemokine receptors in stem cell circulation. Front Biosci 2008; 13: 6820–6833.

31.

Zhang

Chee

Shi

, et al. Intervertebral disc cells produce interleukins found in patients with back pain. Am J Phys Med Rehabil 2016; 95(6): 407–415.

32.

Pattappa

Peroglio

Sakai

, et al. CCL5/RANTES is a key chemoattractant released by degenerative intervertebral discs in organ culture. Eur Cell Mater 2014; 27: 124–136, discussion 36.

33.

Grad

Bow

Karppinen

, et al. Systemic blood plasma CCL5 and CXCL6: potential biomarkers for human lumbar disc degeneration. Eur Cell Mater 2016; 31: 1–10.

34.

Zhang

Nie

Wang

, et al. CCL20 secretion from the nucleus pulposus improves the recruitment of CCR6-expressing Th17 cells to degenerated IVD tissues. PloS One 2013; 8(6): e66286.

35.

Zhang

Liu

Wang

, et al. Production of CCL20 on nucleus pulposus cells recruits IL-17-producing cells to degenerated IVD tissues in rat models. J Mol Histol 2016; 47(1): 81–89.

36.

Sadowska

Touli

Hitzl

, et al. Inflammaging in cervical and lumbar degenerated intervertebral discs: analysis of proinflammatory cytokine and TRP channel expression. Eur Spine J 2018; 27(3): 564–577.

37.

Navone

Marfia

Giannoni

, et al. Inflammatory mediators and signalling pathways controlling intervertebral disc degeneration. Histol Histopathol 2017; 32(6): 523–542.

38.

Molinos

Almeida

Caldeira

, et al. Inflammation in intervertebral disc degeneration and regeneration. J R Soc Interf 2015; 12(104): 20141191.

39.

Gabr

Jing

Helbling

, et al. Interleukin-17 synergizes with IFNγ or TNFα to promote inflammatory mediator release and intercellular adhesion molecule-1 (ICAM-1) expression in human intervertebral disc cells. J Orthop Res 2011; 29(1): 1–7.

40.

Huang

Lin

Chen

, et al. IL-20 may contribute to the pathogenesis of human intervertebral disc herniation. Spine 2008; 33(19): 2034–2040.

41.

Smith

Chiaro

Nerurkar

, et al. Nucleus pulposus cells synthesize a functional extracellular matrix and respond to inflammatory cytokine challenge following long-term agarose culture. Eur Cell Mater 2011; 22: 291–301.

42.

Hoyland

Le Maitre

Freemont

. Investigation of the role of IL-1 and TNF in matrix degradation in the intervertebral disc. Rheumatology 2008; 47(6): 809–814.

43.

Séguin

Pilliar

Roughley

, et al. Tumor necrosis factor-alpha modulates matrix production and catabolism in nucleus pulposus tissue. Spine 2005; 30(17): 1940–1948.

44.

Séguin

Pilliar

Madri

, et al. TNF-alpha induces MMP2 gelatinase activity and MT1-MMP expression in an in vitro model of nucleus pulposus tissue degeneration. Spine 2008; 33(4): 356–365.

45.

Altun

. Cytokine profile in degenerated painful intervertebral disc: variability with respect to duration of symptoms and type of disease. Spine J 2016; 16(7): 857–861.

46.

Schroeder

Viezens

Schaefer

, et al. Chemokine profile of disc degeneration with acute or chronic pain. J Neurosurg Spine 2013; 18(5): 496–503.

47.

Wang

Cai

Shi

, et al. Aging and age related stresses: a senescence mechanism of intervertebral disc degeneration. Osteoarthr and Cartil 2016; 24(3): 398–408.

48.

Kepler

Ponnappan

Tannoury

, et al. The molecular basis of intervertebral disc degeneration. Spine J 2013; 13(3): 318–330.

49.

Hartman

Yurube

, et al. Expression and regulation of metalloproteinases and their inhibitors in intervertebral disc aging and degeneration. Spine J 2013; 13(3): 331–341.

50.

Yang

Zhang

, et al. Intervertebral disc ageing and degeneration: the antiapoptotic effect of oestrogen. Age Res Rev 2020; 57: 100978.

51.

Song

, et al. Estrogen receptor expression in lumbar intervertebral disc of the elderly: gender- and degeneration degree-related variations. Jt Bone Spine 2014; 81(3): 250–253.

52.

Wei

Shen

Williams

, et al. Expression and functional roles of estrogen receptor GPR30 in human intervertebral disc. J Steroid Biochem Mol Biol 2016; 158: 46–55.