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

Background: Inammatory processes exacerbated by IL-1β are believed to be key mediators of disc degeneration and low back pain. However, the underlying mechanism remains unclear. We performed a bioinformatics analysis to identify the key genes that were differentially expressed between degenerative intervertebral disc cells with and without exposure to interleukin (IL)-1β, and explore the related signaling pathways and interaction networks. Methods: The microarray data were downloaded from the Gene Expression Omnibus (GSE 27494). 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 (DAVID), Metascape, Gene Set Enrichment Analysis (GSEA), Search Tool for the Retrieval of Interacting Genes (STRING), Cytoscape, the Venn method, and packages of the R computing language. Results: A total of 260 DEGs were identied, including 161 upregulated genes and 99 down-regulated genes. Gene Ontology (GO) 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 (KEGG) signaling pathway analysis showed that these DEGs were mainly involved in the interactions of cytokine-cytokine receptor interaction, rheumatoid arthritis, tumor necrosis factor (TNF) 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 intervertebral disc degeneration. Conclusions: Bioinformatic analysis showed that CXCL8 and other 9 key genes may play a role in the development of disc degeneration induced by inammatory reactions, and can be used to identify the potential therapeutic target genes.


Materials And Methods
Collection of microarray data The gene pro le GSE27494, which was downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), was produced by the [U133_X3P] Affymetrix Human X3P Array (Platform GPL1352). The gene pro le 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 cells were grown in a 3D collagen construct with or without 10 2 pM IL-1 for a total of 14 days. Following homogenization in TRIzol reagent, the total RNA was isolated and analyzed using microarrays.
Identi cation of differentially expressed genes (DEGs) The limma R package emerged as one of the most widely used statistical tests for identifying DEGs (18). The package was obtained from the 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 drawn using R software.
Functional Annotation of DEGs using a database for annotation, visualization and integrated discovery (DAVID), and Metascape analysis We applied DAVID (https://david.ncifcrf.gov/home.jsp) (version 6.8) (19) to perform Gene Ontology (GO) functional analysis (20) and performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of DEGs (21). The enriched GO terms and pathways with a P-value < 0.05 were identi ed. Moreover, GO and KEGG pathway analyses were performed using Metascape (http://metascape.org/gp/index.html#/main/step1) (22). 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 (GSEA) GSEA primarily analyzes microarray data, using genomic and genetic sequencing to detect the signi cant biological differences in microarray datasets (23). In the present study, we performed GESA 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 as reference gene sets (software.broadinstitute.org/gsea/) (24). Gene set permutations were analyzed 1,000 times. In addition, the normalized enrichment score, normalized P-value, and false discovery rate (FDR) q values were applied to lter the correlative GO terms and pathways.
Construction of the protein-protein interaction (PPI) network, signi cant module, and hub gene network First, we used Metascape (http://metascape.org/gp/index.html#/main/step1) (22) to construct the PPI network and screen the signi cant 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 signi cant. Highly interconnected (hub) genes in the PPI network were analyzed using Cytoscape software (version 3.6.0). The Molecular Complex Detection tool (MCODE) (version 1.5.1) (25) could screen and identify the most signi cant module in the PPI network. Genes with the top 10 highest degrees in the PPI network were viewed as hub genes.

Identi cation of signi cant genes
We delineated a Venn diagram to identify the signi cant common genes among "Metascape_MCODE", "Cytoscape_MCODE", and "Cytoscape_cytoHubba" using Fun Rich software (http://www.funrich.org) (25). Summaries of the function of signi cant genes were obtained via DAVID (https://david.ncifcrf.gov/home.jsp). The R language was used to perform the clustering analysis of signi cant genes based on the gene expression levels.

Screening for DEGs
After the analysis of the datasets (GSE27494) using the limma package, the DEGs between IL-1β-exposed and control disc cells could be presented in a volcano plot (Fig. 1). After setting the threshold value, a total of 260 DEGs were reserved, including 161 upregulated DEGs (Table 1) and 99 down-regulated DEGs (Table 2). The cluster heatmap of these DEGs is shown in Fig. 2.  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 ( Fig. 3A-C). In the KEGG signaling pathway analysis, the 260 DEGs were found to be mainly involved in the interactions of cytokine-cytokine receptor interaction, rheumatoid arthritis, TNF signaling pathway, salmonella infection, chemokine signaling pathway, and other signal transduction pathways (Fig. 3D).
Furthermore, the functional enrichment analysis performed using Metascape showed that the DEGs were signi cantly enriched in NABA MATRISOME ASSOCIATED, extracellular matrix organization, NABA CORE MATRISOME, positive regulation of cell motility, and regulation of cell adhesion ( Fig. 4A-C).
IL-8 is a chemotactic factor that attracts neutrophils, basophils, and T-cells, but not monocytes. It is also involved in neutrophil activation.

CXCL1
C-X-C motif chemokine ligand 1 Has chemotactic activity for neutrophils. May play a role in in ammation and exerts its effects on endothelial cells in an autocrine fashion.
3 CCL20 C-C motif chemokine ligand 20 Chemotactic factor that attracts lymphocytes and, slightly, neutrophils, but not monocytes. Inhibits proliferation of myeloid progenitors in colony formation assays. May be involved in formation and function of the mucosal lymphoid tissues by attracting lymphocytes and dendritic cells towards epithelial cells.

CXCL2
C-X-C motif chemokine ligand 2 Produced by activated monocytes and neutrophils and expressed at sites of in ammation. Hematoregulatory chemokine, which, in vitro, suppresses hematopoietic progenitor cell proliferation.

CXCL5
C-X-C motif chemokine ligand 5 Involved in neutrophil activation.

CXCL3
C-X-C motif chemokine ligand 3 Ligand for CXCR2 (By similarity). Has chemotactic activity for neutrophils. May play a role in in ammation and exert its effects on endothelial cells in an autocrine fashion.

CXCL6
C-X-C motif chemokine ligand 6 Chemotactic for neutrophil granulocytes. Receptor for estrogen.

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 in ammation as the primary etiological factor impacting disc degeneration, and IL-1β was selected as a representative in ammatory cytokine (28). Therefore, we obtained gene expression pro ling data from GEO microarray datasets containing 4 samples of degenerative intervertebral disc cells stimulated with 10 2 pM IL-1β and 4 controls. The key genes that were affected during disc degeneration caused by in ammatory 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. Subsequent analysis con rmed their involvement in important IVDD-related pathways. Our results suggest that these hub genes might be potential biological targets in IVD diagnosis and for the development of therapeutic drugs.
In terms of GO enrichment analysis, we found that most of 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 similarly found that most of DEGs were primarily enriched in the cytokine-cytokine receptor interaction, TNF signaling pathway, chemokine signaling pathway. In addition, CXCL5 and CXCL8 are the most signi cantly upregulated genes in our analysis. In ammation has been correlated with IDD but its role remains controversial. Previous studies have shown that pathological in ammation of IVD and peridiscal space is characterized by increased levels of pro-in ammatory cytokines and chemokines (29)(30)(31). Cytokines and chemokines have three modes of action: (1) recruiting in ammatory cells and activating phagocytosis, (2) stimulating the production of other in ammatory mediators and MMPs, and (3) enhancing matrix degradation (32)(33)(34)(35)(36)(37). And these effects can activate a cascade of detrimental self-promoting events that exacerbate the degeneration of IVD. Furthermore, experimental studies have reported the contribution of pro-in ammatory cytokines in degenerated painful IVD (38, 39).
In addition, our study also showed that the DEGs were mostly enriched in the extracellular region and extracellular matrix. A normal adult IVD is composed of a vascular tissue, which contain a large amount of extracellular matrix (ECM), and the disc cells are responsible for maintaining the integrity of the ECM (40). Changes in the metabolic balance of IVD cells affect the quality and quantity of the ECM and its functional properties, and these changes can, therefore, be related to disc degeneration (41). The underlying molecular mechanisms of IVDD have been previously investigated, and ECM degradation and in ammation have been con rmed to play a critical role in accelerating IVDD progression. Le et al. reported that IL-1β induced the expression of MMPs (16). And the study of Vo et al. found that upregulation of MMP and ADAMTS expression and enzymatic activity is implicated in disc ECM destruction, leading to the development of IDD (42). Our bioinformatics analysis also showed that MMP-3 signi cantly increased in the degenerative intervertebral disc cells stimulated with 10 2 pM IL-1β than controls.
Furthermore, we constructed a PPI network the top 10 hub genes -CXCL8, CXCL1, CCL20, CXCL2, CXCL5, CXCL3, CXCL6, C3, PF4, and GPER1 -were identi ed. It is known that leukocytes, such as neutrophils, macrophages, and lymphocytes, mainly participate in the in ammatory defense response of the body. Chemokines produced by mammalian cells during in ammation are a class of chemotactic and inducible small molecule peptides that are ubiquitous and play an important role in acute and chronic in ammation (43). Based on key cysteine residues involved in disul de bonds, chemokines are classi ed as C, CC, CXC, and CX3C (44). Among the top 10 hub genes, 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.
CXCL8 is a chemotactic factor that attracts neutrophils, basophils, and T-cells. It is involved in neutrophil activation and modulates both acute and chronic in ammation. Ahn et al. found that CXCL8 mRNA expression was associated with the development of radicular pain by back extension and suggested that CXCL8 participates in the pathomechanism of nerve root in ammation in lumbar disc herniations (45). Pederson et al. found that patients with lumbar radicular pain due to disc herniation have increased serum levels of CXCL8 (46). Wang et al. also found higher CXCL8 concentrations in the serum of lumbar disc herniation (LDH) patients (47). A recent study reported that CXCL8 was elevated in the cerebrospinal uid (CSF) of chronic low back pain (LBP) patients with IVDD, compared to pain-free subjects with or without IVDD, and they supported that the IL-8 signaling pathway is a viable therapy for chronic LBP (48). Palada et al. considered that neuroin ammation mediated by elevated CXCL8 concentrations in the CSF and CXCL8 mediated periphery-to-CNS (central nervous system) in ammatory cross-talk contributes to pain in LDH patients (49). Burke et al. reported that the increased level of CXCL8 within the nucleus pulposus may be related to neovascularization in patients with discogenic pain (50). Importantly, the protein expression of CXCL8 is signi cantly increased, concordant with the histological degenerative tissue changes in human NP (51). Similarly, Zhang et al. found that annulus brosus samples from LBP patients had an elevated IL-8 expression compared to controls with scoliosis (52). These experimental results are consistent with the results of our analysis. In the present study, the expression of CXCL8 in IL-1βexposed samples was increased by 7.88-fold compared to control samples, indicating that increased levels of IL-1β affected the expression of the low-back-pain-related factors in intervertebral disc cells, which might become a new target for treating related diseases.
Chemokine CXCL5, which was increased by 9.83-fold in the present study, also called epithelial neutrophil-activating peptide 78 (ENA-78), belongs to the CXC family of chemokines that carry a glutamate-leucine-arginine (ELR) motif and binds to the G-protein CXCR2 to recruit neutrophils and promote angiogenesis (53,54). CXCL5 has been shown to facilitate nociceptive input transmission in the pathogenesis of in ammatory pain (55). Dawes et al. reported that the chemokine CXCL5 is a peripheral mediator of ultraviolet B (UVB)-induced in ammatory pain, likely in humans and rats (56). Xu et al. reported that the upregulation of spinal CXCL5 and CXCR2 is involved in neuropathic pain after nerve injury by regulating GSK-3β activity in a rat model of chronic constriction injury (CCI) of the sciatic nerves (57). In addition, CXCL5 overexpression has been observed in several malignancies, such as pancreatic ductal adenocarcinoma, demonstrating its role in tumor carcinogenesis (58). This nding showed that the role of CXCL5 in intervertebral discs merits further research.
Chemokine CXCL6, also known as granulocyte chemotactic protein 2 (GCP2), recruits in ammatory cells to the site of in ammation by binding to receptors CXCR1 and CXCR2. Sandell et al. reported that the expression of CXCL6 was signi cantly increased in IL-1β-treated human chondrocytes (59). CXCL6 expression in IL-1β-exposed samples was increased by 3.2-fold, which was consistent with the results of a previous study. CXCL6 has also been detected in conditioned medium of induced degenerative discs in organ culture and may contribute to the chemotactic response of induced-degenerative discs (60). Grad  CC chemokine ligand 20 (CCL20), which is expressed by endothelial cells in several tissues, is also known as macrophage in ammatory protein 3α (MIP-3α); it is a 70-amino-acid chemokine that binds exclusively to the chemokine receptor 6 (CCR6) and recruits CCR6expressing cells (68, 69). Zhang et al. reported that NP cells from the extruded and herniated patient group could produce abundant amounts of CCL20 (70). Subsequently, Zhang et al. suggested a potential mechanism for the recruitment of IL-17-producing cells to degenerated intervertebral discs via a CCL20/CCR6 system in vivo, while IL-17 is involved in the auto-immune process of intervertebral disc degeneration in rat models (71).
G protein-coupled estrogen receptor 1 (GPER1) is the receptor of estrogen. There are four majors naturally occurring oestrogens identi ed in women: oestrone (E1), oestradiol (E2), oestriol (E3), and estetrol (E4). In addition, oestrogen receptors consist of two classic nuclear receptors, estrogen receptor (ER)-α and ER-β, and the membrane-bound G-protein-coupled receptor 30 (GPR30) (72). Over the past decade, many researchers have widely studied the relevance between oestrogen and IVDD and found that oestrogen can effectively alleviate IVDD development by inhibiting the apoptosis of IVD cells. Oestrogen can decrease IVD cell apoptosis in multiple ways, including the inhibition of in ammatory cytokines IL-1β and TNF-α, reducing catabolism because of matrix metalloproteinases inhibition, upregulating integrin α2β1 and IVD anabolism, activating the PI3K/Akt pathway, decreasing the oxidative damage, and promoting autophagy (72). Song et al. found that both cytoplasmic and nuclear staining of ERα and ERβ immunoreactivity were observed in nucleus pulposus cells, and ERα and ERβ expression signi cantly decreased, along with the aggravation of IVD degeneration in both males and females (73).
Wei et al. also showed that GPR30 with a high a nity for estrogen is expressed in the human disc NP, can mediate E2 enhanced cell proliferation and in uence disc cell survival (74). In our study, GPER1 was downregulated by 2.42-fold. Therefore, IL-1β decreases the expression of GPER1 in human disc NP, which accelerates disc degeneration.
Chemokine (C-X-C motif) ligand 2 (CXCL2), also known as macrophage in ammatory protein 2 (MIP-2), belongs to the CXC chemokine family, along with growth-regulated protein β and growth-regulated oncogene-2. CXCL2 is 90% identical regarding the amino acid sequence to the related chemokine CXCL1. CXCL2, which is produced in an injured sciatic nerve by partial ligation and is secreted by monocytes, can rapidly recruit neutrophils to sites of in ammation or injured tissues through blood ow, in response to infection or injury (43,75,76). CXCL3 is a small cytokine belonging to the CXC chemokine family, and is also known as the GRO3 oncogene, GRO protein gamma, and macrophage in ammatory protein 2β. PF4, also known as CXCL4, which belongs to the CXCL chemokine family, is predominantly produced by megakaryocytes and α-granules of platelets and has an important role in hemostasis/thrombosis (77). Complement C3 is a central molecule in the complement system and plays a very important role in in ammatory and immune responses. Interestingly, the complement system has been implicated in cartilage degradation, and C3 has been found to be aberrantly increased in the synovial uids from individuals with osteoarthritis and in animal models of osteoarthritis (78, 79). However, we did not retrieve the literature to identify whether IVDD is correlated with CXCL2, CXCL3, C3, and PF4.

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

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.

Funding
Our study was supported by the project of medical discipline priority of Shijingshan District, Beijing.
Authors' contributions NF and SY: design of the study and writing of the manuscript. YH, PD, JL, XK, WZ, YL and LZ: design of the study and revising manuscript critically for important intellectual content. All authors have read and approved the nal manuscript.

Figure 1
Volcano plot presenting the DEGs between IL-1β-exposed and control disc cells. The down-regulated DEGs are marked in green, and the up-regulated DEGs are marked in red.

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 up-regulation and the green color indicates a down-regulation.