Oxidative Stress in Intervertebral Disc Degeneration: New Insights from Bioinformatic Strategies

Oxidative stress has been proved to play important roles in the development of intervertebral disc degeneration (IDD); however, the underlying mechanism remains obscure to date. The aim of this study was to elucidate the vital roles of oxidative stress-related genes in the development of IDD using strict bioinformatic algorithms. The microarray data relevant to the IDD was downloaded from Gene Expression Omnibus database for further analysis. A series of bioinformatic strategies were used to determine the oxidative stress-related and IDD-related genes (OSIDDRGs), perform the function enrichment analysis and protein-protein interaction analysis, construct the lncRNA-miRNA-mRNA regulatory network, and investigate the potential relationship of oxidative stress to immunity abnormality and autophagy in IDD. We observed a significantly different status of oxidative stress between normal intervertebral disc tissues and IDD tissues. A total of 72 OSIDDRGs were screened out for the further function enrichment analysis, and 10 hub OSIDDRGs were selected to construct the lncRNA-miRNA-mRNA regulatory network. There was a very close association of oxidative stress with immunity abnormality and autophagy in IDD. Taken together, our findings can provide new insights into the mechanism research of oxidative stress in the development of IDD and offer new potential targets for the treatment strategies.


Background
It is estimated that approximately 84% adults will suffer from the low back pain (LBP) during their lifetimes [1,2]. Intervertebral disc (IVD) degeneration (IDD) is the principal contributor to the LBP, but the precise pathogenesis remains obscure [3]. Several mechanisms have been proved to be involved in the initiation or progression of IDD, such as mechanical stress, immunity abnormality, metabolic disorders, and oxidative stress [3][4][5][6][7][8]. Especially, oxidative stress has gotten increased attention in the pathogenesis of IDD recently [9][10][11].
Oxidative stress refers to an imbalance between oxidants and antioxidants in the cells and tissues, which will result in the accumulation of reactive oxygen species [6]. Recently, increasing evidence shows that oxidative stress may play vital roles in IDD by inducing the premature senescence and promoting a catabolic phenotype in human nucleus pulposus (NP) cells [9][10][11][12][13][14]. Hu et al. found that exosomes derived from bone mesenchymal stem cells could relive the compression-induced apoptosis of NP cells through inhibiting the oxidative stress [10]. Seol et al. reported that the reduction of oxidative stress by amobarbital was an effective way to prevent the IDD progression [13]. Despite a growing number of relevant studies, the underlying mechanism of oxidative stress in IDD has not been fully elucidated up to now.
Long noncoding RNA (lncRNA) refers to the noncoding RNA whose length exceeds 200 nucleotides, and it can interact with the DNA, RNA, and protein [15]. LncRNA is involved in many diseases via functioning as a competing endogenous RNA (ceRNA) to suppress the microRNA (miRNA) functions and upregulate the targeted gene expression [16][17][18]. However, limited research has been performed to elucidate the ceRNA network of lncRNA in the regulation of oxidative stress in IDD.
Immunity abnormality is a very crucial contributor to the development of IDD [3,19]. The proinflammatory cytokines secreted by the NP cells and the infiltrating immune cells (IICs) are involved in the progression of IDD [5]. Autophagy is an important intracellular degradation process to maintain the homeostasis of intracellular environment by removal of damaged organelles and nonfunction proteins and cyclic utilization of degraded components [20]. The activation of autophagy may be an important protection mechanism to relieve the damage of oxidative stress and thereby prevent the progression of IDD [21,22]. However, the underlying association of oxidative stress with immunity abnormality and autophagy in the development of IDD has not been fully understood.
In this study, we used a series of mature and recognized bioinformatic strategies to determine the following items: (1) The oxidative stress-related genes (OSRGs) involved in the pathogenesis of IDD, their main biological functions, and the potential ceRNA regulatory network; and (2) the possible relationship of oxidative stress to immunity abnormality and autophagy in IDD.

Data Collection and
Processing. This research has been approved by the institutional research ethic committee of Peking University Third Hospital, and the informed consent was exempted because all data was obtained from the common database. The flow chart of this research was shown in Figure 1. Gene expression data of mRNA and miRNA was downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/ geo/). The dataset used for this research was selected according to following criterions: (a) species: Homo sapiens; (b) sequencing methods: microarray data or nextgeneration sequencing data; and (c) samples: the nucleus pulposus in normal IVD tissues and IDD tissues. The following dataset was excluded from this research: expression data from other tissues (e.g., blood sample) and IVD tissues from subjects with scoliosis. Finally, the mRNA expression dataset GSE56081 and miRNA expression dataset GSE63492 were selected and downloaded for further analysis. The expression files of these two datasets were from 5 IDD tissues and 5 normal IVD tissues by the microarray, and the details of these tissue samples were listed in the Supplementary  Table 2). The immunity-related gene list was downloaded from the immunology database and Analysis Portal (ImmPort) database (https://www.immport.org/ home) [24] (Supplementary Table 3

Determination of Oxidative Stress in IDD.
To investigate the potential role of oxidative stress in the pathogenesis of IDD, the oxidative stress score for each tissue sample was calculated using the single sample gene set enrichment analysis (ssGSEA), which was a computational approach to explore whether a priori defined set of genes has statistical significance and concordant differences in two biological conditions for a single sample [26][27][28].
2.3. Identification of Oxidative Stress-Related and IDD-Related Genes (OSIDDRGs) and IDD-Related miRNAs (IDDRmiRNAs). IDD-related genes were obtained using the R package limma package, with the criterion of adjusted p < 0:05 and fold change > 2. The OSIDDRGs were obtained with the intersection of the IDD-related genes and OSRGs using the Venn diagram. IDDRmiRNAs were obtained using the R package limma package, with the standard of p < 0:05 and fold change > 1:20. The volcano plot and heat map were plotted using the R package ggplot2.

Functional Enrichment Analysis of OSIDDRGs. Gene
Ontology (GO) analysis of OSIDDRGs was conducted to interpret their biological process, cell component, and molecular function. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was performed to explore the related signaling pathways of OSIDDRGs. The GO and KEGG analyses were both conducted using the DAVID database (https://david.ncifcrf.gov/) [29], and specific items of GO and KEGG analyses with p <0.05 were selected for further visualization using the R package ggplot2.

Protein-Protein Interaction (PPI) Analysis of OSIDDRGs.
The STRING database (https://cn.string-db.org/) was a widely used tool to explore the association networks of functional proteins [30]. PPI analysis was conducted using the STRING database by imputing the OSIDDRGs into the multiple protein section, and protein pairs with score > 0:40 were selected to construct the PPI network by Cytoscape software (https://cytoscape.org/). The PPI score was calculated using the MCC method by the cytoHubba plug-in, and top 10 OSIDDRGs ranked by the PPI score were selected for further analysis.
2.6. Construction of ceRNA Regulatory Network of lncRNA-miRNA-mRNA. For each OSIDDRG, the correlation analysis between IDDRmiRNAs and hub OSIDDRG was performed using the Pearson test, and significantly related IDDRmiR-NAs with r < −0:7 and p < 0:05 were selected for further intersection with predicted miRNAs. The predicted miRNAs for 10 hub OSIDDRGs were obtained using TargetScanHuman database (http://www.targetscan.org/vert_80/) [31]. The intersection of significantly related IDDRmiRNAs and predicted miRNAs was performed using the Venn diagram to obtain the candidate miRNAs for the construction of miRNA-mRNA pairs. Then, the upstream targeted lncRNAs for candidate miRNAs were predicted to construct the lncRNA-miRNA pairs using the Predicted model in the LncBase Predicted database (http://carolina .imis.athenainnovation.gr/diana_tools/web/index.php?r= lncbasev2/index-predicted) [32]. Ultimately, the ceRNA 2 Oxidative Medicine and Cellular Longevity regulatory network of lncRNA-miRNA-mRNA was constructed and presented by the Cytoscape software.

Correlation Analysis between Oxidative
Stress and Immunity Abnormality. The immunity score for each tissue sample was calculated using the ssGSEA algorithm based on the immunity-related gene list [5] and compared between normal IVD tissues and IDD tissues. The immune infiltration analysis was performed using the CIBERSORT database (http://cibersort.stanford.edu/) [33]. The difference of IICs was compared between normal IVD tissues and IDD tissues, and IICs with p < 0:05 were considered as IDD-related IICs (IDDRIICs). The correlation analysis between oxidative stress score and immunity score and the correlation analysis between 10 hub OSIDDRGs and 8 types of IDDRIICs were conducted to explore the potential relationship between oxidative stress and immunity abnormality in IDD using the Pearson test. To further elucidate the interaction between oxidative stress and immunity abnormality in the development of IDD, the miRNA-mRNA-IDDRIIC regulatory network was built by integrating the IDDRmiRNA-OSIDDRG pairs with r < −0:70 and p < 0:05 with OSIDDRG-IDDRIIC pairs with |r | >0:70 and p < 0:05. Ultimately, the miRNA-mRNA-IDDRIIC regulatory network was constructed and presented using the Sankey diagram.

Correlation Analysis between Oxidative Stress and
Autophagy. The autophagy score for each tissue sample was calculated using the ssGSEA algorithm based on the autophagy-related gene list [26] and compared between normal IVD tissues and IDD tissues. To obtain the autophagy-related and IDD-related genes (ATGIDDRGs), the intersection between autophagy-related genes and IDD-related genes was conducted using the Venn diagram. Then, the PPI analysis of ATGIDDRGs was conducted to obtain the top 10 hub ATGIDDRGs using the MCC method by Cytoscape software. Furthermore, the correlation analysis between hub OSIDDRGs and hub ATGIDDRGs was conducted.
2.9. Statistical Analysis. Any analysis in this study was conducted by the R software 4.1.2. The ssGSEA score and percentage of IICs between normal IVD tissues and IDD tissues were compared using the Student's t test and p < 0:05 indicating there was a statistical difference between two groups. Correlation analysis in this study was performed using the Pearson test, and the correlated pair with |r | >0:70 and p < 0:05 was considered as statistically significant.

Confirmation of Differential Oxidative Stress
Score between Normal IVD Tissues and IDD Tissues. The oxidative stress score was calculated by ssGSEA algorithm based on the expression of OSRGs, and there was a significant difference of oxidative stress score between normal IVD tissues and IDD tissues (p < 0:01) (Figure 2(a)). PCA analysis showed the expression of OSRGs could clearly distinguish the normal IVD tissues and IDD tissues (Figure 2(b)). Therefore, oxidative stress might play vital roles in the pathophysiology of IDD.   Table 5).

Correlation Analysis between Oxidative Stress and
Autophagy in IDD. An obvious difference of autophagy score was observed between the two groups (p = 0:049) (Figure 7(a)), which manifested autophagy was involved in the development of IDD. There was a significant relationship between oxidative stress score and autophagy score (r = 0:80, p < 0:01) (Figure 7(b)), which indicated there might be a potential link between oxidative stress and autophagy in the pathogenesis of IDD. Eighteen ATGIDDRGs were obtained with the intersection between IDD-related genes and autophagy-related genes (Figure 7(c)), and top 10 hub ATGIDDRGs were selected using PPI analysis (Figure 7(d)). Then, the correlation analysis between hub OSIDDRGs and hub ATGIDDRGs was performed and 74 OSIDDRG-ATGIDDRG pairs were observed (Figure 7(e)). The TXNRD1-HSPA5 pair was the most positively correlated pair (r = 0:99, p < 0:01) (Figure 7(f)), and MAPK1-CDKN1A pair was the most negatively correlated pair (r = −0:92, p < 0:01) (Figure 7(g)).

Discussion
IDD has become a leading cause for LBP, which results in huge economic burden for the patients and medical care

18
Oxidative Medicine and Cellular Longevity system [1][2][3]. Despite abundant researches in recent years, the definite etiology for IDD remains obscure to date [3,5,22]. Recently, accumulating evidence showed that oxidative stress might play important roles in the pathogenesis of IDD, but the underlying mechanism was still unclear [12,13,34]. In this study, for the first time, we determined the different status of oxidative stress between normal IVD tissues and IDD tissues using the ssGESA algorithm, which indicated oxidative stress had a very vital role in the IDD. Furthermore, we screened out 72 OSIDDRGs, identified their relevant biological functions, and constructed the ceRNA regulatory network using a series of bioinformatic strategies. Moreover, we also preliminarily elucidated the association of oxidative stress with immunity abnormity and autophagy in IDD. Our findings provided new insights into the researches about the oxidative stress in the development of IDD. Through several strict bioinformatic algorithms, 10 hub OSIDDRGs were identified, including IL6, PRDX1, MCL1, HMOX1, TXNRD1, MAPK1, HIF1A, FOXO1, JUN, and JAK2. The IL6 was a famous inflammatory cytokine and could promote the initiation and progression of IDD [35]. The HMOX1 could alleviate the senescence of NP cells through inducing the autophagy [36]. The MAPK1 has been reported to be upregulated by EZH2 and promote the IDD via suppressing the miR-129-5p [37]. The HIF1A could alleviate the compression-induced apoptosis of NP-derived stem cells trough upregulating the autophagy via the HIF1A-BNIP3-ATG7 axis [38]. The FOXO1 was a vital regulator of IVD homeostasis through the direct regulation of autophagy, adaptation to hypoxia, and resistance to oxidative stress [39]. The c-Jun was found to improve the expression of TGF-β and promote the NP cell proliferation via reducing the apoptosis and inflammatory response [40]. The JAK2/STAT signaling path-way has been testified to protect the IVD from oxidative stress induced degeneration [41]. To our knowledge, the biological roles of PRDX1, MCL1, and TXNRD1 in IDD have not been investigated up to now. The PRDX1 has been verified to be involved in the regulation of NK cells and protective autophagy in hepatocellular carcinoma [42]. The MCL1 was a famous antiapoptotic protein and could stimulate the progression and drug resistance of thymic carcinoma [43]. The TXNRD1 was a vital ferroptosis-related gene and has been approved to promote the invasion, progression, and metastasis of hepatocellular carcinoma [44]. Based on previous studies [42][43][44], we made the assumption that PRDX1, MCL1, and TXNRD1 might participate in IDD progression via the regulation of oxidative stress, cell apoptosis, and immunity status, but future studies should be conducted to interpret the precise underlying mechanism. Moreover, to further explore the biological functions of OSIDDRGs, we performed the function enrichment analysis and identified several crucial KEGG pathways which might be involved in IDD. In fact, several of these identified signaling pathways have been verified in IDD, including TNF signaling pathway [45], estrogen signaling pathway [7], MAPK signaling pathway [46], HIF-1 signaling pathway [47], Toll-like receptor signaling pathway [48], and glutathione metabolism [49]. However, how the OSIDDRGs participate in IDD via these signaling pathways has not been fully understood and needs the further investigation.
LncRNA can affect the gene expression by competitively binding miRNAs, namely, ceRNA mechanism, which has been proved to play important and complicated roles in the initiation or progression of human diseases, including IDD [15,50,51]. Lan  the potential molecular mechanisms of miRNAs, lncRNAs, and circular RNAs in the progression of IDD were described based on the latest literatures, which suggested that noncoding RNAs could serve as potential targets for the IDD treatment [53]. Moreover, a bioinformatic analysis conducted by Wang et al. showed that miRNAs could function as novel targets for preventing and treating IDD by regulating their target genes [54]. In a similar manner, a meta-analysis conducted by Sherafatian et al. also observed the vital potential roles of abnormal miRNAs in the development of IDD [55]. In the current study, we constructed a strict ceRNA regulatory network of lncRNA-miRNA-mRNA, which was composed of 63 lncRNAs, 9 candidate miRNAs, and 6 hub OSIDDRGs, such as TMEM191C-miRNA-1184-FOXO1 axis and C1orf143-miRNA-4747-5p-TXNRD1. These ceRNA pathways shed new lights on the investigation of oxidative stress in IDD, which are worth our validating and investigating in the future. Immunity abnormality is another crucial contributor to the development of IDD, and many inflammatory cytokines and IICs have been identified to aggravate the IDD [3,5,19]. In our study, we confirmed the differential immunity status between normal IVD tissues and IDD tissues using the ssGSEA algorithm. Subsequently, we identified 8 types of IDDRIICs and performed the correlation analysis among them. We discovered that neutrophils-B_cell_memory pair and dendritic_cells_activated-NK_cells_activated pair were the most negatively and positively correlated pair, respectively. These significantly correlated pairs were worth well of further investigation because few studies have focused on the interaction among IICs in IDD up to now. To further elucidate the potential mechanism of the link between oxidative stress and immunity abnormality in IDD, like the ceRNA regulatory network of lncRNAs, we constructed the miRNA-mRNA-IDDRIICs regulatory network, such as miR-1184-FOXO1-neutrophils axis and miR-4533-MCL1-B_ cells_memory axis. We hypothesized these miRNA-mRNA pairs might serve vital roles by changing the percentage of IDDRIICs in the pathogenesis of IDD, which deserved the special attention in future researches.
Previous studies have indicated that oxidative stress might promote the IDD via the regulation of autophagy [22,34,56]. Tang et al. used the hydrogen-peroxide to induce the oxidative stress in IDD, and they observed an obviously increased expression of LC3-II protein, indicating that oxidative stress could promote the formation of autophagosome in NP cells [34]. Similarly, Xu et al. study showed that a treatment of 400 μM tert-butyl-hydroperoxide for 6 h could also induce autophagy in NP cells by increasing the expression of LC3-II and Beclin-1 [57]. In our study, there was a significant correlation between oxidative stress score and autophagy score, suggesting a potential interaction between oxidative stress and autophagy in the pathomechanism of IDD. We identified 10 hub ATGIDDRGs through the PPI analysis and then attempted to interpret their relationship with OSIDDRGs using the Pearson analysis. Interestingly, most of hub OSIDDRGs had a significant association with hub ATGIDDRGs, indicating a close relationship between oxidative stress and autophagy in IDD. For instance, MAPK1-CDKNA1 pair was the most negatively correlated pair, which has not been fully investigated to date and deserved our major attention in the future.
Several limitations should be considered when interpreting our results. First, different from tumors, only few sequencing files of IDD were uploaded into public databases, which might be attributed to the rareness of control normal IVD tissues. As a result, the sample size of this study was relatively small, which might affect the stringency of our findings. Second, although the bioinformatic strategies used were strict and widely recognized, all findings in this study were obtained only based on the microarray data from GEO database, which needed further experimental verification. Third, although we tried to elucidate the association of oxidative stress with immunity abnormality or autophagy in IDD, the analysis performed in the current study was relatively thin, and in-depth scientific researches were needed. Forth, the retrospective design of this research would inevitably induce the bias, which might affect the reliability of results; therefore, prospective researches with the strict design were necessary to confirm our findings in the future.
Taken together, we found that normal IVD tissues and IDD tissues could be distinctly distinguished via the expression of OSRGs, and there was a significantly different oxidative stress score between the two groups, which indicated that oxidative stress might play important roles in the development of IDD. Then, through a series of strict bioinformatic strategies, the biological functions of OSIDDRGs were elucidated and 10 hub OSIDDRGs were identified (e.g., IL6, PRDX1, and MCL1). Moreover, the ceRNA network of lncRNA-miRNA-mRNA was constructed to further interpret the potential mechanism of oxidative stress in IDD. Furthermore, the relationship of oxidative stress to immunity abnormality and autophagy in the development of IDD was preliminarily evaluated. Our findings could provide new insights into the oxidative stress in IDD and deliver potential targets for the treatment of IDD.