ABSTRACT
ABSTRACT
Background
Bioinformatics analysis was performed on gene expression profile microarray data to identify the key genes activated through the TNF-α/TNFR1 signaling pathway in intervertebral disc degeneration (IDD). The common differentially expressed genes (co-DEGs) were calculated in nucleus pulposus (NP) cells and annulus fibrosus (AF) cells under TNF-α treatment or TNFR1 knockdown, which reveals the potential mechanism of TNF-α involvement in IDD and may provide new therapeutic targets for IDD.
Methods
Differentially expressed genes (DEGs) in TNF-α-treated or TNFR1-knockdown NP cells and AF cells were identified. Further analysis of the gene ontology (GO), signaling pathways and interaction networks of the DEGs or co-DEGs were conducted using the Database for Annotation, Visualization and Integrated Discovery, STRING Database, and Cytoscape software. The relationship between genes and musculoskeletal diseases, including IDD, was assessed with the Comparative Toxicogenomics Database. The predicted microRNAs corresponding to the co-DEGs were also identified by microRNA Data Integration Portal.
Results
In NP cells, the DEGs (|log2FoldChange|>2, adj.P < 0.01) were identified including 48 DEGs by TNF-α treatment and 74 DEGs by TNFR1 knockdown; in AF cells, correspondingly, 105 DEGs were identified. The co-DEGs between NP cells and AF cells were calculated including CXCL8, ICAM1, BIRC3, RELB, NFKBIA, and TNFAIP3. They may be the hub genes that were significantly associated with both NP cells and AF cells through the TNF-α/TNFR1 signaling pathway. The co-DEGs and corresponding predicted miRNAs may be potential therapeutic targets for IDD.
Conclusions
CXCL8, ICAM1, BIRC3, RELB, NFKBIA, and TNFAIP3 may have a synergistic effect on TNF-α-induced IDD development.
Abbreviations: IDD: Intervertebral disc degeneration; NP: Nucleus pulposus; AF: Annulus fibrosus; co-DEG: Common differentially expressed gene; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; PPI: Protein-protein interaction.
Background
Intervertebral disc degeneration has developed into a global social health problem that often causes chronic low back pain and a socioeconomic burden1,2. The etiologies of IDD are multifaceted, including genetic susceptibility, abnormal biomechanical stress, aging, inflammatory factors, and so on1-4. Secreted by NP cells, AF cells, macrophages, T cells and neutrophils, inflammatory factors play an important role in IDD and discogenic pain5,6. Inflammatory factors activate the receptors in NP cells and AF cells, regulate the transcription or translation of downstream genes, and thus, participate in the regulation of the extracellular matrix, which in turn causes a series of pathogenic reactions5.
Among these cytokines, tumor necrosis factor (TNF)-α is one of the most important cytokines in IDD. TNF-α signals through its receptors TNFR1 and TNFR2, thereby affecting downstream genes and participating in the process of IDD7. TNF-α is highly expressed in degenerated intervertebral disc tissue and is deeply involved in multiple pathological processes of IDD, including matrix degradation, inflammatory response, apoptosis, autophagy, and cell proliferation8. TNF-α promotes the expression of matrix-degrading enzyme genes, such as the matrix metalloproteinases MMP-3 and MMP-13, a disintegrin and a metalloproteinase with the thrombospondin motifs ADAMTS4, ADAMTS-5, ADAMTS-9, and ADAMTS-15, together exacerbating IDD9. In addition, it has been reported that anti-TNF-α therapy can alleviate IDD and low back pain to some extent10-12.
TNFR1, also known as TNFRSF1A, encodes a protein that is a member of the TNF receptor superfamily. It is one of the main receptors of TNF-α, which can activate the transcription factor NF-κB and play an essential role in the regulation of inflammation13-15. TNF-α/TNFR1 signaling regulates the expression of downstream genes, of which the upregulated genes were considered pivotal genes to exacerbate inflammation in the intervertebral disc and may be potential biomarkers or therapeutic targets for IDD.
The aim of our study was to identify differentially expressed genes (DEGs) in NP cells with or without TNF-α treatment or TNFR1 knockdown and those in AF cells with or without TNF-α treatment to provide new biological targets for IDD.
Methods
Sample preparation & cell culture
The datasets GSE132958 and GSE41883 were downloaded from the NCBI GEO Datasets, and the expression profiling arrays were generated using GPL16791 Illumina HiSeq 2500 (Homo sapiens) and GPL1352 [U133_X3P] Affymetrix Human X3P Array. In the dataset of GSE132958, the NP tissue was obtained from surgical tissue of patients with DDD-related back pain (n = 5). All patients were accompanied by axial back pain, and confirmed by MRI, their intervertebral discs showed signs of degeneration16. According to the description of GSE132958, the NP tissue was dissected, minced and digested enzymatically. Later, the cells were counted and monolayer- cultured in DMEM with 10% FBS. And then, the NP cells were tested for knockdown of TNFR1 with dCas9-KRAB system before RNA-seq analysis. The TNR1 knockdown NP cells were treated with or without 10 ng/mL of TNF-α for 24 hours. In the dataset of GSE41883, human disc tissue samples were obtained from patients with herniated discs and degenerative disc disease. AF cells were expanded for use in 3D in a collagen sponge with or without 10e3 pM TNF-α for 14 days. To further verify the result by bioinformatic analysis, additional real-time PCR experiments were conducted in NP cells and AF cells. And in the real-time PCR experiments, the normal human NP cells which were purchased from ScienCell, USA, and the AF cells which were obtained from surgical tissue, were cultured in DMEM (HyClone, USA) with 10% FBS (Gibco, USA), and maintained in a humidified cell incubator (Thermo Fisher Scientific, USA) with 5% CO2, 21% O2 and 74% N2 at 37 °C. The cells were treated with TNF-α (50 ng/ml) for 24 h after serum starvation for 6 h.
Data processing
The normalized count files of the GSE132958 dataset were directly downloaded from the NCBI GEO database. The raw data of GSE41883 were obtained, and R packages of “affy” and “affyPLM” were used to assess the dataset. The Benjamini-Hochberg method was used to adjust the original p-values, and the false discovery rate (FDR) was used to calculate the fold-changes (FC). The R package of “limma” was used for further analysis, and the DEGs were screened using R with the following data filtering settings: log2FC>2 or -log2 FC>2 and adjusted P-value<0.01.
In NP cells, the significant DEGs (|log2FoldChange|>2, adj.P < 0.01) in the TNF-α stimulation group (abbreviation: TNF-α group) compared to those in the control group were collected. Then, the significant DEGs in the TNFR1-knockdown group in the condition of TNF-α (abbreviation: TNFR1-knockdown group) in contrast to those in the TNF-α stimulation group alone were collected. Additionally, Venn diagrams were made for the common DEGs between the TNF-α group and the TNFR1 knockdown group. These common DEGs may account for the destructive effect of TNF-α/TNFR1 in human NP cells. Similarly, the significant DEGs in TNF-α-stimulated AF cells were collected. Finally, we calculated and made Venn diagrams for co-DEGs in human NP cells and AF cells under specific conditions.
The R package of ggPlot2 and pheatmap were used for visualization of the DEGs collected.
Identification of protein-protein interaction (PPI) networks of DEGs
The DEG lists were uploaded to the STRING version 11.0 online analysis platform (https://string-db.org), and the PPI networks were obtained for the DEGs. The network and tab-separated values were further analyzed by Cytoscape version 3.7.1 (http://cytoscape.org/) and inserted Apps MCODE (Molecular Complex Detection v1.5.1).
Functional annotation for DEGs
The DEG lists were uploaded to the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 (http://david.ncifcrf.gov/); then, GO, including the cell components, molecular functions, and biological processes, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs were analyzed and visualized by the R package of ggPlot2. The Cytoscape app, ClueGo, was also used to conduct functional annotation of the biological processes of the DEGs in AF cells and visualize the GO terms and KEGG pathways.
Identification of co-DEGs associated with IDD and predicted miRNAs
The Comparative Toxicogenomics Database (http://ctdbase.org/) (CTDbase) was used to further analyze the relationship between final co-DEGs and IDD. A gene Venn diagram was created to show the inferred associations with IDD. MicroRNA Data Integration Portal (http://ophid.utoronto.ca/mirDIP/index.jsp) (MirDIP) was used to analyze and predict miRNAs targeting the co-DEGs.
RNA isolation & real-time PCR
Total RNA was isolated from human NP cells with TRIzol reagent (TAKARA, Dalian, China) according to the manufacturer’s instructions. cDNA was synthesized by reverse transcription with PrimeScript RT reagent Kit (TAKARA, Dalian, China). Real-time PCR was performed using TB Green Premix Ex Taq II (Takara, Dalian, China). Bio-Rad CFX Connect (Bio-Rad, USA) and Bio-Rad CFX Manager were used to perform real-time PCR and analyze the data. The results were quantified using the ΔΔCt method. All primers were synthesized by Generay, Shanghai, China. The primer sequences used for real-time PCR were as followed (Table 1.)
Table 1. Primer sequences used for real-time PCR
Statistical analysis
The real-time PCR data are presented as the mean ± standard deviation (SD) of three experiments. GraphPad Prism 7 (GraphPad Software Inc. La Jolla, CA) was used to perform the statistical analysis. Data were analyzed via the t-test if normally distributed and the Mann-Whitney test or if not normally distributed. Values of p < 0.05 were considered statistically significant.
Results
Identification of DEGs of NP and AF
Compared to the normal control cells, GSE132958 NP cells contained a total of 277 DEGs in the TNF-α group that were identified by screening; these genes included 196 upregulated genes and 81 downregulated genes (|log2FC|>2, adj.P < 0.05) (Figure 1(a)). As shown in Table 2 and Figure 1.C, 48 genes were markedly differentially expressed in the TNF-α group (log2FC>2, adj.P < 0.01) including 29 upregulated and 19 downregulated DEGs, which were used for further analysis and gene annotation. To investigate how TNF-α/TNFR affected the process of disc degeneration, we also identified the DEGs in the TNFR1 knockdown group following TNF-α treatment. Compared to the TNF-α treatment alone group, the knockdown group had a total of 335 DEGs (log2FC>2, adj.P < 0.05) (Figure 1(b)), of which 74 genes were markedly differentially expressed (|log2FC|>2, adj.P < 0.01) including 24 upregulated and 50 downregulated DEGs. As shown in Table 3 and Figure 1(d), 50 genes were markedly downregulated when TNFR1 was knocked down by the lentiviral CRISPR epigenome editing system under the condition of TNF-α. These DEGs were considered to be associated with the mechanism of TNF-α, which explained its destructive effect on the intervertebral disc through the TNF-α/TNFR1 signaling pathway.
Published online:
26 August 2020Figure 1. Hierarchical clustering analysis of TNF-α-related and TNFR1 knockdown-related differentially expressed genes in NP cells. (A) TNF-α-related DEGs (|log2FC|>2, adj.P<0.05) and (B) TNFR1 knockdown-related DEGs were plotted in the volcano plot. The heat maps of (C) the DEGs between TNF-αtreatment and control (log2FC>2, adj.P<0.01) and (D) the DEGs between TNFR1 knockdown + TNF-α treatment and TNF-αtreatment only (log2FC>2, adj.P<0.01) were shown correspondingly
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/icts20/2021/icts20.v062.i05/03008207.2020.1797709/20210908/images/medium/icts_a_1797709_f0001_oc.jpg"})
Figure 1. Hierarchical clustering analysis of TNF-α-related and TNFR1 knockdown-related differentially expressed genes in NP cells. (A) TNF-α-related DEGs (|log2FC|>2, adj.P<0.05) and (B) TNFR1 knockdown-related DEGs were plotted in the volcano plot. The heat maps of (C) the DEGs between TNF-αtreatment and control (log2FC>2, adj.P<0.01) and (D) the DEGs between TNFR1 knockdown + TNF-α treatment and TNF-αtreatment only (log2FC>2, adj.P<0.01) were shown correspondingly
Table 2. Top 30 markedly differentially expressed genes in TNF-α-treated NP cells
Table 3. Top 30 markedly differentially expressed genes in TNFR-knockdown NP cells with TNF-α treatment
To reveal the role of TNF-α in AF cells, further analysis was conducted in TNF-α treated cells. As was reported in GSE41883 dataset, 105 DEGs (|log2FC|>2, adj.P < 0.01) were identified in AF cells17. The common and core DEGs in NP cells and AF cells were finally identified, which may be potential gene biomarkers and therapeutic targets of TNF-α related IDD.
Functional GO terms and pathway enrichment analysis
The GO functional annotation analysis of the 48 DEGs in the TNF-α group showed that inflammatory response (P = 2.60E-07), cellular response to lipopolysaccharide (P = 1.60E-04) and interleukin-1 (P = 6.90E-04) and immune response (P = 5.10E-04) were the top 4 biological altered genes (Figure 2(a)). KEGG pathway enrichment analysis also revealed that these genes mainly participated in the TNF signaling pathway (P = 1.60E-09), NF-kappa B signaling pathway (P = 3.90E-07), rheumatoid arthritis (P = 4.10E-07), NOD-like receptor signaling pathway (P = 1.10E-06), and cytokine-cytokine receptor interaction (P = 1.40E-04) (Figure 2(b)).
Published online:
26 August 2020Figure 2. GO terms and KEGG pathway enrichment. (A)(B) GO term enrichment and KEGG pathway of the 48 DEGs relative to TNF-α treatment. (C)(D) GO term enrichment and KEGG pathway of the 74 DEGs relative to TNFR1 knockdown. Dot sizes represent counts of enriched DEGs, and dot colors represent negative log2 P values
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/icts20/2021/icts20.v062.i05/03008207.2020.1797709/20210908/images/medium/icts_a_1797709_f0002_oc.jpg"})
Figure 2. GO terms and KEGG pathway enrichment. (A)(B) GO term enrichment and KEGG pathway of the 48 DEGs relative to TNF-α treatment. (C)(D) GO term enrichment and KEGG pathway of the 74 DEGs relative to TNFR1 knockdown. Dot sizes represent counts of enriched DEGs, and dot colors represent negative log2 P values
Responsively, related to the dysregulated genes in TNFR1 knockdown group with TNF-α treatment in NP cells, the top 4 terms of biological processes also included inflammatory response (P = 7.00E-12), immune response(P = 4.10E-07), and response to lipopolysaccharide response (1.70E-07) as well as NF-kappaB signaling (P = 7.50E-05) (Figure 2(c)). Similarly, extracellular space (P = 1.50E-12) and chemokine activity (P = 2.70E-05) were also markedly related to cellular components and molecular functions of these dysregulated DEGs in TNFR1 knockdown group. The top 4 terms of KEGG pathway enrichment were identical to those of TNF-α treatment group.
The functional GO terms and pathway enrichment analysis in AF were referred to the previous bioinformatics analysis18.
PPI network analysis
The common DEGs in the TNF-α group and e TNFR1-knockdown group were identified and calculated as shown in the Venn diagram (Figure 3(a)). The PPI network of these 24 common DEGs was then created as shown in Figure 3(b) and further analyzed by MCODE. The hub nodes included BIRC3, CCL2, RELB, TNFAIP3, ICAM1, CXCL8, NFKBIA and MAP3K14, which were considered to be the hub genes of TNF-α that induced a destructive effect on NP cells (Figure 3(c)). Based on a previous study, 105 DEGs in TNF-α-stimulated AF cells were also identified from the GSE41883 dataset. The PPI network revealed the potential relationships among these 105DEGs in AF cells (Figure 3(d)). Similarly, MCODE was used to further predict the inner relationship among 16 DEGs shown in Figure 3(e).
Published online:
26 August 2020Figure 3. Venn diagrams and PPI network for the DEGs in NP and AF respectively. (A) Venn diagrams of the DEGs relative to both TNF-α treatment (log2FC>2, adj.P<0.01) and TNFR1 knockdown (log2FC>2, adj.P<0.01) in NP cells. (B) PPI networks from the DEGs in NP cellsconstructed using STRING database. (C) MOCDEsuggested a possible relationship among the DEGs in NP cells. (D) PPI networks from the DEGs relative to TNF-α treatment in AF cells. (E) MCODE suggested a possible relationship among the DEGs in AF cells
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/icts20/2021/icts20.v062.i05/03008207.2020.1797709/20210908/images/medium/icts_a_1797709_f0003_oc.jpg"})
Figure 3. Venn diagrams and PPI network for the DEGs in NP and AF respectively. (A) Venn diagrams of the DEGs relative to both TNF-α treatment (log2FC>2, adj.P<0.01) and TNFR1 knockdown (log2FC>2, adj.P<0.01) in NP cells. (B) PPI networks from the DEGs in NP cellsconstructed using STRING database. (C) MOCDEsuggested a possible relationship among the DEGs in NP cells. (D) PPI networks from the DEGs relative to TNF-α treatment in AF cells. (E) MCODE suggested a possible relationship among the DEGs in AF cells
Co-DEGs of NP cells and AF cells
The DEGs in both NP cells and AF cells were identified and the common DEGs were then calculated. A Venn diagram was created and shown in Figure 4(a). The co-DEGs for NP and AF included CXCL8, BIRC3, RELB, ICAM1, NFKBIA, TNFAIP3, WTAP and SLC39A14. These co-DEGs were then compared to genes with inferred associations to IDD on the Comparative Toxicogenomics Database, further narrowed down to the following gene list: CXCL8, BIRC3, ICAM1, TNFAIP3, RELB and NFKBIA (Figure 4(b)). The PPI network by Cytoscape also verified the close relationship of the key genes (Figure 4(c)). ClueGO was used to visualize the functional annotation of the biological processes and KEGG pathways and it turned out that these co-DEGs were significantly related to NFκB signal pathway (Figure 4(d)). Real-time PCR also revealed that the mRNA expression levels of these 6 co-DEGs were significantly enhanced by TNF-α in human NP cells and AF cells (Figure 4(e,f)).
Published online:
26 August 2020Figure 4. Venn diagrams and PPI network. (A) Venn diagrams of the key genes in NP cells relative to TNF-α/TNFR1 and key genes in AF relative to TNF-α. Co-expressed genes included NFKBIA, TNFAIP3, CXCL8, RELB, ICAM1, BIRC3, SLC39A14 and WTAP. (B) Venn diagrams of the co-expressed genes of NP and AF and IDD-related genes constructed using CTDbase. The inferred genes included NFKBIA, TNFAIP3, CXCL8, RELB, ICAM1 and BIRC3. (C) PPI networks from the inferred genes in (B) constructed using STRING database for DEGs and clueGO confirmed the close relationship among these 6 genes. (D) ClueGo suggested the potential GO and KEGG pathway of these 6 co-DEGs. (E) The mRNA expression levels of NFKBIA, TNFAIP3, CXCL8, RELB, ICAM1 and BIRC3 were significantly increased in TNF-α-stimulated NP cells and (F) AF cells. Data represent the mean ± standard deviation of three independent experiments. *P < 0.05
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/icts20/2021/icts20.v062.i05/03008207.2020.1797709/20210908/images/medium/icts_a_1797709_f0004_oc.jpg"})
Figure 4. Venn diagrams and PPI network. (A) Venn diagrams of the key genes in NP cells relative to TNF-α/TNFR1 and key genes in AF relative to TNF-α. Co-expressed genes included NFKBIA, TNFAIP3, CXCL8, RELB, ICAM1, BIRC3, SLC39A14 and WTAP. (B) Venn diagrams of the co-expressed genes of NP and AF and IDD-related genes constructed using CTDbase. The inferred genes included NFKBIA, TNFAIP3, CXCL8, RELB, ICAM1 and BIRC3. (C) PPI networks from the inferred genes in (B) constructed using STRING database for DEGs and clueGO confirmed the close relationship among these 6 genes. (D) ClueGo suggested the potential GO and KEGG pathway of these 6 co-DEGs. (E) The mRNA expression levels of NFKBIA, TNFAIP3, CXCL8, RELB, ICAM1 and BIRC3 were significantly increased in TNF-α-stimulated NP cells and (F) AF cells. Data represent the mean ± standard deviation of three independent experiments. *P < 0.05
Identification of co-DEGs associated with IDD and predicted miRNAs
The co-DEGs of NP cells and AF cells with or without TNFR1 knockdown or TNF-α treatment were then put into CTDbase to predict or verify the association with musculoskeletal diseases,including IDD. According to the results of CTDbase, the 6 co-DEGs most targeted muscular diseases and arthritis diseases, and they also showed a close relation to IDD (Figure 5). Prediction analysis using mirDIP was used to identify the predicted microRNAs targeting the co-DEGs, as listed in Table 4.
Published online:
26 August 2020Figure 5. Relationship of the co-expressed genes to musculoskeletal diseases including IDD based on the CTD database
![]({"type":"image","src":"/na101/home/literatum/publisher/tandf/journals/content/icts20/2021/icts20.v062.i05/03008207.2020.1797709/20210908/images/medium/icts_a_1797709_f0005_oc.jpg"})
Figure 5. Relationship of the co-expressed genes to musculoskeletal diseases including IDD based on the CTD database
Table 4. Top 5 predicted miRNAs of the co-DEGs of NP and AF
Discussion
Currently, a great number of people suffer from low back pain, which is often related to IDD19. TNF-α, one of the most common inflammatory factors in IDD, exerts its function mainly through TNFR1, and by activating the downstream genes, TNF-α accelerates the process of IDD6,13–15,20-22 .Therefore, the identification of key genes downstream of the TNF-α/TNFR1 pathway in NP cells and AF cells may provide new insights into potential therapeutic targets for IDD.
Our study was conducted to identify the co-DEGs of NP cells and AF cells under the condition of TNF-α or TNFR1 knockdown. The introduction of TNF-α in NP cells resulted in the differential expression of 48 genes while knockdown of TNFR1 led to the differential expression of 74 genes and thus, the common genes accounted for the destructive effect of TNF-α through TNFR1 in NP. The DEGs were markedly enriched in GO terms associated with inflammatory and immune response, response to lipopolysaccharide and extracellular space as well as chemokine activity, which explained and indicated that TNF-α significantly resulted in the pro-inflammatory effect in the extracellular matrix catabolism in NP cells. Liu et al. also identified the DEGs in human AF cells exposed to TNF-α and IL-1β and had similar results with our analysis of GSE41883 that DEGs in AF cells were also markedly enriched in terms associated with inflammatory response and extracellular space23.
6 co-DEGs were finally identified in the analysis including CXCL8, TNFAIP3, NFKBIA, RELB, ICAM1 and BIRC3. GO and KEGG analysis in Figure 4(d) revealed that these DEGs focus on the NF-kappa B signaling pathway while cytoplasmic pattern recognition receptor, nucleotide-bind oligomerization domain containing (NOD), nucleotide-binding domain, leucine rich repeat containing receptor (NLR) and toll-like receptor 4 (TLR4) signaling pathway significantly participated in the pathological process of these core genes. This suggested that activated NF-κB pathway mainly accounted for the destructive effect of TNF-α in IDD. MCODE was then used to enable searches for dense clique-like structures and find clusters (highly interconnected regions) in the network, and a cluster among CXCL8, TNFAIP3, NFKBIA, RELB, ICAM1 and BIRC3 was suggested. It was considered that there could be a synergy among these 6 core genes in the pathological process of TNF-α/TNFR related IDD. However, no study was reported on the relationship among these genes and it was worth of further research.
CXCL8, a member of the CXC chemokine family, also known as IL-8, is a major mediator of the inflammatory response and increases concordant with histological degenerative tissue changes24. The CXCL8 inhibitor also showed a protective effect on chronic radicular neuropathic pain because of disc herniation25. CXCL8 exists as an inflammatory factor and participates in the activation of NF-κB signal pathway, which explains the differential expression of CXCL8 and may contribute to IDD. As for ICAM1, also termed CD54, it was identified as a classical senescence markers of disc cells for the first time and it was upregulated in senescent NP cells under the conditions of low glucose, hypoxia, high osmolality and the absence of serum26. Besides, the addition of IFN-γ, TNF-α or IL-17 also demonstrated a marked increase in ICAM-1 expression both in human NP cells and AF cells27. In addition, it was recently reported that IFN-γ and LPS promoted IDD by CCL2/ICAM1 activation through the FAK/ERK/GSK3 and PKCδ signaling pathways in AF cells28. Above all, ICAM1 was considered to show destructive effect on IDD. These results revealed that knockdown or inhibition of CXCL8 or ICAM1 may demonstrate a protective effect on disc degeneration.
KEGG pathway enrichment analysis revealed that NF-κB pathway was markedly associated with the biological process of these co-DEGs. Moreover, RELB, an NF-kB subunit, participates in the noncanonical NF-κB pathway, transfers into the nucleus with p52 and regulates the expression of targeted genes29. However, NFKBIA encodes a member of the NF-kappa-B inhibitor family and inhibits the activity of dimeric NF-κB/REL complexes by trapping REL dimers in the cytoplasm by masking their nuclear localization signals30. The differential expression of RELB and NFKBIA also verifies the important role of NF-κB pathway in TNF-α/TNFR1 related IDD.
Interestingly, BIRC3, a member of the IAP family of proteins, inhibits apoptosis by binding to TRAF1 and TRAF2, while TNFAIP3, whose expression is rapidly induced by the TNF, could also inhibit TNF-mediated apoptosis31-33. The differential expression of BIRC3 and TNFAIP3 may explain the limited effect of TNF-α on apoptosis induction in the intervertebral disc. Recently it was reported that BIRC3 was identified as a co-DEG while NF-κB signaling as the common affected pathway in the juvenile and adult Kashin-Beck disease (KBD) and the OA34. Some studies also showed that TNFAIP3 could bind to the WD40 domain of ATG16L1 to control the autophagic response and NF-κB activation in the intestine and that it also restricted MTOR signaling and promoted autophagy in CD4 + T cells, which proposed a new role of TNFAIP3 in the regulation of autophagy and inflammation35,36. Since little research has been conducted on these 2 genes in the disc, further analysis is necessary to reveal the role and mechanism of BIRC3 and TNFAIP3 in IDD.
Besides, according to the CTDbase, these 6 co-DEGs had a close relation with several musculoskeletal diseases, including muscular diseases, rheumatoid arthritis, muscle weakness and intervertebral disc degeneration. This suggests the co-DEGs may play a role in the process of these diseases, and the relations with intervertebral disc degeneration deserve further study. Plenty of evidence revealed that microRNAs were recognized as important regulators of gene expression regulating NP cell functions and that they may be a novel strategy of biological therapy for IDD37. In addition, our previous study also showed that microRNA-145 overexpression increased the matrix synthesis of NP cells38. Therefore, mirDIP database was used to identify the corresponding predicted microRNAs for these 6 genes and according to the integrated scores, it was suggested that miR-106b-5p, miR-223-3p, miR-374A-5p, miR-19b-3p, miR-7-5p and miR-381-3p could be potential targets for biological therapy for IDD.
Taken together, PPI network analysis indicated a possible synergy of CXCL8, ICAM1, BIRC3, TNFAIP3, RELB and NFKBIA in the TNF-α/TNFR1-mediated NF-κB pathway, which may participate in the initiation, course, and termination of TNF-α function. Therefore, the function and regulation of each hub gene, as well as the in-depth relationship among these genes, needs further research and analysis. The corresponding predicted miRNAs of these hub genes are also worth investigating to further and fully understand the role and mechanism of TNF-α in the development of IDD.
There were several limitations in this study. First, TNFR1 knockdown by the lentiviral CRISPR epigenome editing system was conducted on NP cells but not AF cells. No relative public data could be found in AF cells and we only integrated NP cell in the condition of TNF-α or TNFR1 knockdown and AF cell in the condition of TNF-α. Second, bioinformatics analysis of DEGs in NP or AF cells exposed to TNF-α or IL-1β had been conducted so we focus on the DEGs of TNF-α/TNFR1 related IDD for the first time. Last, further experimental research is needed to verify the role and mechanism of these co-DEGs in the pathological process of IDD.
Conclusion
The hub genes of CXCL8, ICAM1, BIRC3, TNFAIP3, RELB, and NFKBIA may be significantly associated with IDD. The top 5 miRNAs for each may be potential biomarkers or therapeutic targets for TNF-α-related IDD, especially miR-106b-5p, miR-223-3p, miR-374A-5p, miR-19b-3p, miR-7-5p and miR-381-3p. Thus, this will provide new insight into the mechanism and therapy of IDD.
| Species | Gene | Sense (5ʹ to 3ʹ) | Antisense (5ʹ to 3ʹ) |
|---|---|---|---|
| Homo | ACTB | CTACCTCATGAAGATCCTCACCGA | TTCTCCTTAATGTCACGCACGATT |
| CXCL8 | ACTGAGAGTGATTGAGAGTGGAC | AACCCTCTGCACCCAGTTTTC | |
| RELB | CAGCCTCGTGGGGAAAGAC | GCCCAGGTTGTTAAAACTGTGC | |
| ICAM1 | ATGCCCAGACATCTGTGTCC | GGGGTCTCTATGCCCAACAA | |
| BIRC3 | GCCATTGACTTTTCTGTCGCC | GCAAAGCAAGCCACTCTGT | |
| TNFAIP3 | TCCTCAGGCTTTGTATTTGAGC | TGTGTATCGGTGCATGGTTTTA | |
| NFKBIA | CTCCGAGACTTTCGAGGAAATAC | GCCATTGTAGTTGGTAGCCTTCA |
| Gene | logFC | P-Value | Adj. P-Value | Gene | logFC | P-Value | Adj. P-Value |
|---|---|---|---|---|---|---|---|
| CACNG4 | −3.339392578 | 1.32E-08 | 2.61E-05 | GSAP | 2.430991799 | 1.89E-05 | 1.70E-03 |
| RNASE7 | 2.876176602 | 1.87E-08 | 2.61E-05 | CSF1 | 2.571999903 | 2.83E-05 | 2.42E-03 |
| LOC101927653 | −2.665526619 | 2.24E-08 | 2.61E-05 | C15orf48 | 6.898390681 | 2.89E-05 | 2.45E-03 |
| CBWD3 | −2.399189555 | 2.86E-08 | 2.61E-05 | RELB | 2.87187701 | 2.98E-05 | 2.52E-03 |
| LOC102724467 | −2.267010759 | 3.27E-08 | 2.61E-05 | GDF6 | −2.00961386 | 3.25E-05 | 2.73E-03 |
| LINC00629 | −2.267010759 | 3.27E-08 | 2.61E-05 | IFNAR2 | 2.614254909 | 3.56E-05 | 2.97E-03 |
| NFKBIA | 2.648153841 | 1.79E-06 | 1.90E-04 | PPP1R1B | −3.172042977 | 4.05E-05 | 3.31E-03 |
| SLC39A14 | 2.391039758 | 2.25E-06 | 2.36E-04 | TNIP1 | 2.074165643 | 4.36E-05 | 3.53E-03 |
| SOD2 | 3.452904962 | 4.12E-06 | 4.14E-04 | AKR1B1 | 2.962601894 | 4.49E-05 | 3.62E-03 |
| ICAM1 | 4.953132858 | 6.48E-06 | 6.43E-04 | MAL2 | −2.164615134 | 4.79E-05 | 3.81E-03 |
| TNFRSF9 | 4.067323206 | 8.00E-06 | 7.81E-04 | MAP3K14 | −2.018734225 | 5.00E-05 | 3.97E-03 |
| HR | −2.256917133 | 8.09E-06 | 7.81E-04 | RRAD | 2.943712121 | 5.42E-05 | 4.23E-03 |
| IL6 | 3.682230007 | 1.20E-05 | 1.12E-03 | PITPNM3 | −2.40352538 | 5.57E-05 | 4.30E-03 |
| ADAMTS15 | −5.456018122 | 1.24E-05 | 1.16E-03 | WTAP | 2.138745644 | 5.74E-05 | 4.41E-03 |
| CCL2 | 3.563990463 | 1.55E-05 | 1.42E-03 | CXCL8 | 7.871770356 | 6.03E-05 | 4.61E-03 |
| Gene | logFC | P-Value | Adj. P-Value | Gene | logFC | P-Value | Adj. P-Value |
|---|---|---|---|---|---|---|---|
| NRXN2 | −3.243361569 | 9.55E-09 | 2.09E-05 | SLC39A14 | −2.269334398 | 4.10E-06 | 5.30E-04 |
| PTGIR | −2.275700529 | 2.19E-08 | 2.09E-05 | LGR5 | 2.189204979 | 5.59E-06 | 7.02E-04 |
| CABP1 | −2.241476566 | 2.26E-08 | 2.09E-05 | MSX2 | −3.901020302 | 7.48E-06 | 9.21E-04 |
| DLEU1-AS1 | 2.241476566 | 2.26E-08 | 2.09E-05 | NFKB2 | −2.350621833 | 9.07E-06 | 1.11E-03 |
| GPR132 | 2.241476566 | 2.26E-08 | 2.09E-05 | LCE1 C | −5.505302205 | 1.18E-05 | 1.41E-03 |
| PSG7 | 2.241476566 | 2.26E-08 | 2.09E-05 | MFAP2 | 2.324777206 | 1.27E-05 | 1.51E-03 |
| HIST1H2AI | −2.241476566 | 2.26E-08 | 2.09E-05 | KCCAT211 | −2.88897809 | 1.45E-05 | 1.69E-03 |
| PCOLCE-AS1 | 2.241476566 | 2.26E-08 | 2.09E-05 | IFNAR2 | −2.314466003 | 1.55E-05 | 1.78E-03 |
| LOC101927815 | 2.241476566 | 2.26E-08 | 2.09E-05 | CCL2 | −3.891527927 | 1.62E-05 | 1.83E-03 |
| LUZP2 | 4.414101147 | 4.91E-07 | 6.92E-05 | PITPNM3 | 2.525457087 | 1.66E-05 | 1.86E-03 |
| GDF6 | 2.09585347 | 5.21E-07 | 7.31E-05 | NLGN4Y | −2.244129081 | 1.90E-05 | 2.08E-03 |
| INHBB | 2.564248835 | 7.23E-07 | 1.00E-04 | BIRC3 | −3.959180461 | 1.93E-05 | 2.12E-03 |
| HR | 2.564769296 | 1.61E-06 | 2.16E-04 | TNFAIP3 | −3.956311857 | 2.06E-05 | 2.23E-03 |
| IL10RB-AS1 | −3.693969412 | 3.75E-06 | 4.90E-04 | C1QTNF1 | −2.30040955 | 2.13E-05 | 2.27E-03 |
| CSF1 | −2.436624758 | 4.06E-06 | 5.27E-04 | EFHD1 | 3.267270028 | 2.16E-05 | 2.29E-03 |
| Gene | Predicted miRNAs | Gene | Predicted miRNAs |
|---|---|---|---|
| CXCL8 | hsa-miR-106b-5p | RELB | hsa-miR-7-5p |
| hsa-miR-106a-5p | hsa-miR-204-5p | ||
| hsa-miR-17-5p | hsa-miR-18a-5p | ||
| hsa-miR-20b-5p | hsa-miR-211-5p | ||
| hsa-miR-20a-5p | hsa-miR-491-5p | ||
| TNFAIP3 | hsa-miR-19b-3p | ICAM1 | hsa-miR-223-3p |
| hsa-miR-19a-3p | hsa-miR-940 | ||
| hsa-miR-23b-3p | hsa-miR-299-3p | ||
| hsa-miR-18b-5p | hsa-miR-483-3p | ||
| hsa-miR-23a-3p | hsa-miR-431-5p | ||
| NFKBIA | hsa-miR-381-3p | BIRC3 | hsa-miR-374a-5p |
| hsa-miR-300 | hsa-miR-369-3p | ||
| hsa-miR-449a | hsa-miR-374b-5p | ||
| hsa-miR-449b-5p | hsa-miR-142-3p | ||
| hsa-miR-548 c-3p | hsa-miR-155-5p |
Acknowledgments
We appreciated Jianchao Sun for helping with the data processing of the manuscript.
Disclosure statement
No potential conflict of interest was reported by the author(s).
Availability of data and materials
The datasets used during the current study are available in the NCBI GEO Datasets, STRING, DAVID, CTDbase and MirDIP databases..
Declarations
The authors report no conflicts of interest. The authors alone are responsible for the content and writing of this manuscript.
Ethics approval and consent to participate
This study was approved by the ethics committee of Sun Yat-Sen Memorial Hospital, Sun Yat-Sen University, and by the Animal Care and Use Committee of Sun Yat-Sen University.