DOI: https://doi.org/10.21203/rs.3.rs-2096080/v1
As a natural process of aging, intervertebral disc degeneration is more prone to degeneration, with limited repair ability, which is closely related to ageing and excessive manual labor. The main manifestations include the formation of fissures in the annulus fibrosus, the reduction of water in the intervertebral disc, and the decrease of elasticity. However, the molecular mechanism of intervertebral disc degeneration is still unclear.
In this study, key biomarkers in intervertebral disc degeneration were identified through bioinformatics. GSE70362 and GSE56081 were downloaded from the GEO database, and then the two datasets were differentially analyzed and validated for their expression, immune infiltration analysis, functional enrichment analysis, and potential drug prediction through the Connectivity Map (CMap) database.
A total of 352 and 9815 differential genes were identified by GSE70362, GSE56081, respectively. The up-regulated and down-regulated genes of the two datasets were intersected with ferroptosis genes to obtain five key genes that were significantly correlated with immune cell content, namely AKR1C3, CKB, KRT19, MT1G and MUC1. The ROC results showed that the five core genes could well predict the occurrence and development of the disease. In addition, the results of CMap suggested that four drugs, including 1-Phenylbiguanide, LY-2183240, Flubendazole and Penciclovir, have the potential to reverse intervertebral disc degeneration.
Exploring the expression levels of five key genes in intervertebral disc degeneration is conductive to providing new ideas for the prevention and treatment of intervertebral disc degeneration. Moreover, Flubendazole and Penciclovir have the potential to provide options for clinical treatment of intervertebral disc degeneration.
The symptom of low back pain is extremely common among all age groups [1]. Low back pain has a high disability rate, and it also causes high economic burden [2; 3; 4]. It is reported that the economic loss caused by poor work performance of patients is even higher than medical expenditure[5]. The degeneration of the intervertebral discs (IVDs) is a major contributor to low back pain. Despite the high mortality and disability rate, effective drug treatment and intervention are still lacking.
Previous studies reported the increased levels of inflammatory cytokines, increased aggrecan and collagen degradation in IVD cell phenotypes[6]. IVD is acknowledged to be mediated by the abnormal production of proinflammatory cytokines and catabolic molecules secreted by both nucleus pulposus (NP) cells and annulus fibrosus (AF) cells, as well as macrophages, T cells and other leukocytes[7; 8; 9; 10]. Immune-related genes and immune cells infiltration act as an essential role in IVD development and progression. Therefore, the comprehensive analysis of the immune infiltration may provide a new reference for the treatment and prognosis of IVD.
Ferroptosis is a type of regulated necrosis which is driven by the accumulation of lipid peroxidation and lethal reactive oxygen species[11]. In this study, the role of ferroptosis-related genes in IVD was comprehensively evaluated by integrating the two GEO databases. Besides, biological characteristics and immune landscape were compared between the low- and high-expression gene groups through a broad array of bioinformatic analyses. In conclusion, deeper insights into the molecular changes and pathogenesis of disc degeneration throughout the process were successfully provided, including immune infiltration and prediction of potential therapeutic drugs, and the findings provided a reference for immune cells that may drive the progression of disc degeneration, which contributed to the understanding of the occurrence and progression of IVD.
Gene Expression Data Acquisition
The GSE63061 Series Matrix File (Platform: GPL10295) was downloaded from the NCBI Gene Expression Omnibus (GEO) public database (http://www.ncbi.nlm.nih.gov/geo/). The GSE63061 dataset contains 24 groups of patients’ expression profiles, including eight healthy controls and 16 IVD patients. The GSE85426 Series Matrix File (Platform: GPL15314) was obtained from the GEO public database, involving 10 groups of patients’ gene expression profile data with five cases of normal group and five cases of IVD patients. Totally, 428 ferroptosis-related genes were collected from the GeneCards database (https://www.genecards.org/).
Functional Enrichment Analysis
Function annotation of DEGs was performed using the Metascape database (www.metascape.org) to comprehensively explore the functional associations of them. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway of specific genes were carried out. The differences were statistically significant at Min overlap≥3 and p≤0.01.
Analysis of Immune Cell Infiltration
The infiltrating scores of 16 immune cells and the activities of 13 immune-related pathways were further evaluated with the single-sample gene set enrichment analysis (ssGSEA) in the "gsva" R package. Risk score and immune cell infiltration were calculated using Spearman correlation analysis. P < 0.05 was considered to be statistically significant.
Enrichment analysis for transcription factors
Transcription initiation eukaryotes, a highly complex process, often requires the assistance of various protein factors. Transcription factors and RNA polymerase II form a transcription initiation complex and participate in the process of transcription initiation together. According to their functions, transcription factors can be divided into two categories. The first category is general transcription factors. When they form a transcription initiation complex with RNA polymerase II, transcription can start at the correct position. A cis-acting element is a sequence flanking a gene that affects gene expression. Cis-acting elements include promoters, enhancers, regulatory sequences, and inducible elements, and they participate in the regulation of gene expression. The cis-acting element itself does not encode any protein, but only provides important binding site to interact with the trans-acting factor. This analysis was conducted using “cisTarget” R package [12]. In this process, “mm9-500bp-upstream-7species.mc9nr.feather” was applied in the Gene-motif rankings database.
Gene Set Enrichment Analysis (GSEA)
Gene set enrichment analysis (GSEA) is a computational method that uses priori gene sets and ranks of genes to identify significantly changed biological themes in microarray datasets [13]. In order to explore the molecular mechanism of core genes, GSEA was performed to compare the differences in signaling pathways between the high-expression group and the low-expression group. Gene set permutations were executed 1000 times for each analysis. The criteria of significantly enriched pathways were normalized (p value <0.05).
Statistical Analysis
All statistical analyses and plotting were performed using R (Version 4.0). The differences were statistically significant at the level of P < 0.05 for two-sided tests.
Identification of Potential Small Molecule Drugs
The Connectivity Map (CMap) is an online resource potentially capable of discovering new drugs and stimulating new hypotheses of connection among diseases, genes, and therapeutic drugs. Gene expression profile data processed by 1309 small drug molecules were collected from the CMap. The treatment conditions are varied, including different drugs, different concentrations, different treatment durations, and so on. In this study, the differentially expressed genes of the disease were used to predict the targeted therapy drugs for the disease.
Identification of DEGs
The GSE70362 dataset was downloaded from the NCBI GEO public database, and the total number of patients was 24 (IVD group, 16; control group, 8). Differentially expressed genes between the IVD and control groups in the GSE70362 datasets were identified using the Limma R package. Genes retained from diverse tissues were selected using the following criteria: a | log2 (fold-change) | > 0.585 and adjusted P-value < 0.05. A total of 352 DEGs, including 173 downregulated and 179 upregulated, genes were identified (Figure 1a).
Functional enrichment of DEGs
To further explore the biological functions and signaling pathways associated with the 352 DEGs, DEGs were analyzed through GO and KEGG. The significant GO terms and KEGG pathway are shown in Figure 2. As shown in the GO analysis, the DEGs were remarkably enriched in “skeletal system development”, “regulation of cell adhesion”, “innate immune response”, “complement and coagulation cascades” and “extracellular structure organization”.
Identification of ferroptosis DEGs and PPI analysis
The dataset, including 428 ferroptosis-related, was also obtained from GeneCards database and intersected with GSE70362 to identify ferroptosis DEGs. A total of nine up-regulated and eight down-regulated genes were found altogether. The Venn diagram of the DEGs is shown in Figure 3a. STRING was used to construct ferroptosis DEGs PPI networks with confidence scores ≥0.4 and visualized using Cytoscape software. Totally, 16 nodes and 30 edges were identified (Figure 2).
Verification of Key Genes Associated with Ferroptosis
The GSE56081 dataset was downloaded from the NCBI GEO public database. Then, the expression trend of these 17 ferroptosis genes in GSE56081 was analyzed. Five (one up-regulated gene; four down-regulated genes) ferroptosis genes were found to have expression trend consistent with the GSE70362 dataset (Figure 1).
Immune Infiltration Analyses
The microenvironment mainly incorporates immune cells, extracellular matrix, a variety of growth factors, inflammatory factors, and special physicochemical properties, which significantly affects the diagnosis of diseases and the sensitivity of clinical treatment. By investigating the association between key genes and immune infiltration in disease data sets, the potential molecular mechanisms that influences the progression of IVD were further explored. The results indicated that the proportion for B cells and Th1 cells in the IVD group were remarkably lower compared with those in the control groups, while the proportion of many components were significantly higher than those in the control groups, such as Parainflammation, Th2 cells and Type II IFN Response. The interaction between the immune cells is shown in Figure 3.
Analysis of key Genes and Immune Infiltration
The relationship between key genes and immune cells was further explored. The results revealed that multiple key genes were highly correlated with immune cells. The AKR1C3 was positively correlated with Type II IFN Response, but negatively related to T cell co inhibition, B cells and APC co stimulation. Besides, CKB was positively correlated with APC co inhibition, T cell co inhibition, and Type I IFN Response, but negatively correlated with Type II IFN Response. KRT19 had a positive relationship with APC co inhibition. MT1G was positively correlated with Type II IFN Response and tumor-infiltrating lymphocytes (TIL). MUC was positively correlated with T cell co inhibition and B cells, but negatively correlated with Type II IFN Response (Figure 4). Next, the correlations between the five key genes and different immune factors, including chemokines, immunomodulator and receptors, from the TISIDB database were acquired (Figure 5). The results of these analyses proved the key genes were closely associated with immune cell infiltration and played a significant role in the immune microenvironment.
The Association Between Key Genes and Disease-related Genes
A total of 964 disease-related genes were obtained from the GeneCards database (https://www.genecards.org/). Top 20 disease-related genes were identified based on the relevance score from the Genecards database. The relevance of five key genes and top 20 disease-related genes in expression levels were calculated. It was found that the expression levels of the five key genes were significantly correlated with several disease-related genes. (Figure 5).
Enrichment analysis for transcription factors
Motif enrichment analysis was conducted using the five key genes, and it was discovered that they are regulated by the mechanisms related to several transcription factors. Transcription factors enrichment analysis was performed using Area Under Curve (AUC). Motif-TF and selection of important genes were also carried out. The analysis results proved that the Motif with the highest normalized enrichment score (NES: 6.65) was annotated as cisbp_M1838. A total of three genes, namely CKB, KRT19 and MUC1, were enriched in this motif. All enriched motifs and corresponding transcription factors of key genes have been displayed in Figure 6.
Functional Analysis of the Key Genes
Then, the specific signaling pathways behind the five key genes were further investigated. GSEA was implemented to explore the potential molecular mechanisms, thus leading to the progression of IVD. AKR1C3 was enriched for “positive regulation of cyclase activity”, “positive regulation of translation” in GO analysis and “glycan biosynthesis and glycan degradation” in KEGG analysis. CKB was enriched for “peptidyl asparagine modification and “costamere” in GO analysis and “linoleic acid metabolism” in KEGG analysis. KRT19 was enriched for “histone H3 K27 methylation” and “cytosolic ribosome” in GO analysis and “ECM receptor interaction” and “ribosome” in KEGG analysis. MT1G was enriched for “iron ion transmembrane transport” and “protein membrane adaptor activity” in GO analysis and “adherens junction” and “arachidonic acid metabolism” in KEGG analysis. MUC1 was enriched for “myelin maintenance” and “acyl-CoA binding” in GO analysis and allograft rejection and apoptosis in KEGG analysis (Figure 7).
Related Drugs Screening for IVD Treatment
Subsequently, 57 miRNAs and 106 mRNA-miRNA relationship pairs from the five core genes were acquired through the miRcode database and the results were visualized using Cytoscape (Figure 6). Moreover, the CMap database was used to predict the target drug of the differentially expressed genes form the GSE70362 dataset. The results illustrated that the expression profiles of 1-Phenylbiguanide, LY-2183240, Flubendazole, and Penciclovir perturbation-induced gene expression had the most significantly negative correlation with the expression profiles of disease perturbations, which indicated that drugs could alleviate or even reverse the course of disease (Figure 8).
In the present study, the differences in gene expression between IVD patients and healthy controls were examined, and it was found that the changes in immune and metabolic functions are associated with these genes. AKR1C3, CKB, KRT19, MT1G and MUC1 were identified to have a therapeutic and preventive role in IVD, which may be a marker for early diagnosis and prognosis.
The heterogeneity of immune microenvironment in IVD is very high, and the type and number of infiltrating immune cells vary greatly in different locations. In this study, ssGSEA algorithms results revealed a higher percentage of Th2 cells and a lower percentage of Th1 cells in IVD groups. According to this outcome, it was found that the Th1/Th2 ratio of IVD patients were lower than that of the healthy control. In normal circumstances, Th1 and Th2 cytokines maintain a relative balance, which is essential for maintaining an immune equilibrium. An imbalance of Th1/Th2 cell may interfere with the cytokine network and lead to the occurrence and development of many diseases [14; 15].
IFN II response is significantly higher in the IVD group in our study, which is consistent with previous studies [10; 16; 17]. IFN-γ is a critical regulatory protein for both innate and adaptive immunity and plays an important role in immunological cell signaling [18]. Macrophages induce a large number of inflammatory responses during autoimmune and autoinflammatory diseases. TNF, IL-1α, IL-1β, IL-6, IL-8, IL-17, IL-2 and IL-10 are associated proinflammation mediators in IVD degeneration [10; 19; 20; 21]. Almost all of these inflammatory responses are exacerbated by interferon-gamma (IFN-γ), which increases the macrophage production of inflammatory mediators [22]. It has been found that the immune microenvironment containing macrophages or macrophage phenotypes involved in incomplete healing of annulus fibrosus and accelerated the progression of IVD [23]. Gabr, M. A. [24] et al. reported that IL-17 and IFN-γ contribute to the release of inflammatory mediators, and an increase in intercellular adhesion molecule 1 (ICAM − 1) expression. Meanwhile, it was discovered that the DEGs was correlated to multiple pathways, including cell adhesion, innate immune response, ECM-receptor interaction and adherens junction.
Through the verification of GSE56081 dataset, five key genes were finally identified: AKR1C3, CKB, KRT19, MT1G and MUC1. As a member of aldo-keto reductase family, AKR1C3 has a capacity to regulate the cell proliferation and differentiation in a hormone-independent manner [25]. A dual function for AKR1C3 in the metabolism of exogenous compounds was found. The aberrant expression of AKR1C3 is associated with the progression and aggressiveness of various diseases, including polycystic ovary syndrome, obesity, hyperandrogenemia, and tumors [26; 27; 28; 29]. In prostate cancer, it was confirmed that AKR1C2 and AKR1C3 mediates prostate cell proliferation via prostaglandin conversion through PI3K/Akt signaling pathways [30]. Previous studies revealed that PI3K/AKT signaling pathway can regulate the major extracellular component, aggrecan in NP cells and increase the activity of intervertebral disc cell matrix synthesis [31; 32; 33]. These findings indicated that AKR1C2 may contribute to the development of IVD via PI3K /AKT pathway. Additionally, our results displayed an inverse association between the expression of AKR1C3 and COL9A3. COL9A3 gene encodes one of the three alpha chains of type IX collagen that plays a connective role in creating cross-links between annulus fibrosis region and nucleus pulposus [34; 35; 36]. Paassilta et al. reported that at least one Trp3 allele was at three-fold greater risk of developing IVD [37]. Toktas et al. [38] found that patients with COL9A3 Trp3 were more likely to obtain a higher Pfirrmann score. Further studies should be conducted to investigate the potential pathway in AKR1C3 regulating IVD.
Creatine kinase B (CKB), a brain-type creatine kinase, reversibly catalyzes the transfer of phosphate from phosphocreatine to ADP to generate creatine and ATP [39]. CKB is crucial for the tissue with high fluctuating energy demands and plays an important part in distribution and local supply of ATP [40]. Downregulated expression of CKB severely affects cell energy homeostasis, thus leading to different cancers, obesity, and other metabolic diseases [41; 42; 43; 44]. In the model of CKB-knockout mice, Chang et al. observed markedly reduced bone resorption activity in osteoclasts [45]. The finding indicated the bone loss effect of CKB in osteoclasts [45]. What’s more, CKB was positively associated with Vitamin D receptor (VDR) gene that is a member of the nuclear steroid nuclear hormone receptor superfamily of ligand-inducible transcription factors [46]. The vitamin D endocrine system is essential to the calcium resorption, normal bone mineralization and remodeling [46]. The disorder of vitamin D metabolism is associated with pathological conditions of articular cartilage and intervertebral disc tissue, in particular osteoarthritis (OA) [47; 48] and IVD [49; 50; 51]. Furthermore, VDR is one of the most widely studied candidate genes that are supposed to involve IVD. Polymorphisms in VDR gene has been demonstrated to be associated with IVD in Japanese and Finnish populations. Sansoni et al. illustrated that FokI polymorphisms of the VDR gene, particularly the FF genotype and F allele, presented approximately 2-fold risk factors in the process of disc herniation [52]. Besides, L. Chen et al. performed a meta-analysis and found that VDR FokI polymorphism of VDR may be correlated to IDD among Caucasians [53]. These findings may provide an experimental basis for the treatment of IVD.
KRT19 belongs to a family of keratins, which is responsible for structural rigidity and multipurpose scaffolds [54]. KRT19 is often used as a positive marker for the identification of nucleus pulposus cells [55]. Moreover, some studies have shown that KRT19 expression in degenerated NP tissue is lower than that in normal NP tissue [55; 56; 57], which is consistent with our outcome that the expression of KRT19 is lower in the IVD tissues. Lv et al. [56] found that KRT19 presented a high NP:AF ratio and a high NP:AC ratio in degenerated NP. However, little is known about the potential mechanism about the downregulation of KRT 19 in aging IVD tissues or cartilage tissue.
Metallothioneins (MTs) are a family of low-molecular weight and cysteine-rich intracellular proteins that are expressed in response to oxidative stress and metal-induced toxicity [58]. MTs might also regulate the function of zinc-dependent proteins, such as p53, transcription factor, and chemo-sensitivity, by donating or taking away of zinc ions. MT1G is a widely studied isoform of MTs in humans, and it has been reported to play an important part in repressing carcinogenesis and inhibiting tumor metastasis and cancer cells differentiation because of excessive oxidative damage. The expression of MT1G is closely associated with aggressive phenotype, drug resistance, and poor prognosis for tumor progression of different types [59; 60; 61; 62].
It has been well established that MUC1, a membrane-bound mucin expressed on the surface of airway epithelial cells [63; 64], also serves an anti-inflammatory role in the course of various infections [65; 66; 67]. Some studies have demonstrated that the suppression of MUC1 expression can induce a recruitment of macrophages and trigger a hallmark of overstimulated immune responses [67; 68; 69]. McAuley et al. [68] found that the MUC1-knockout (KO) mice exhibited more susceptible to influenza A virus. They observed that MUC1-KO mice have a higher number of macrophages and a high level of pro-inflammatory cytokines, particularly IL-6 and MCP-1. Moreover, under the repetitive airway Pseudomonas aeruginosa (Pa) infection, MUC1-KO mice presented significantly elevated macrophages and airspace enlargement compared with Muc1 wild-type mice [70]. Based on our results, the course of IVD is closely associated with immune response and MUC1 may have vital mechanism behind IVD. Gruber et al. [71] ever reported that the expression of MUC1 increased in high-grade IVD. Interestingly, in vitro studies showed a significant downregulation of MUC1 in human annulus cells while they were exposed to IL-1 and TNF, which indicated that MUC1 may participate important immune response and result in the consumption of it. On the other hand, it was speculated that degenerated intervertebral may have different stages and we have reasons to believe that the activation of immune response may exist in some specific stages of IVD. Promoting the expression of MUC1 in IVD may serve as a new direction in the treatment of IVD.
Since there were no effective drugs for IVD, the online database was widely used to aid our understanding of human diseases and the prediction of new drugs. CMap is an easily accessible and effective tool in exploring new drugs, and it has been confirmed by many studies. From the CMap database, four compounds, including 1-Phenylbiguanide, LY2183240, Flubendazole and penciclovir, were identified and may have significant potential therapeutic effects on IVD. Penciclovir is a novel acyclic nucleoside analogue that has been shown to be effective against HSV1 and HSV2 [72]. In addition, penciclovir can alleviate neuroinflammation in neurodegenerative diseases [73], and inflammation can promote the degradation of extracellular matrix of intervertebral disc, thereby leading to intervertebral disc dysfunction and structural damage [74]. The role of penciclovir in IVD degeneration deserves further discussion. Flubendazole was approved in 1980 for the treatment of gastrointestinal nematode infection [75]. Recent studies have reported that flubendazole has significant antiproliferative activity in vitro and in vivo [76]. Inducing autophagy cell death seems to be the key to flubendazole mediated TNBC cell growth inhibition [76], and the aging and degeneration of IVDs are related to autophagy signal pathway transduction[77].
Some limitations of our study should be mentioned. First, the performance of our five key genes should be evaluated in more IVD datasets. Second, the main results of our study were based on public datasets, so they require further validation by actual experiments. Third, the clinical information of the IVD sample is unknown, which may influence the results of our study.
In conclusion, the five immune-related genes signature were identified to be related to IVD, with potentially substantial clinical significance. Most importantly, two prospective drugs for the treatment of IVD were identified. More studies are needed to confirm the therapeutic targets of these genes. Furthermore, bioinformatic analyses will provide new approaches for studying the pathogenesis and treatment of IVD degeneration.
Data Available
Data utilized for supporting the results of this study are available on the GEO database (https://www.ncbi.nlm.nih.gov/geo/).
Author contribution
Jian-guo Fang and Duo-jun Wang designed the study. Jian-guo Fang and Cai Liu wrote the manuscript and analyzed the data. Zai-Jun Lin: supervision, reviewing and editing. All authors have read and agreed to the published version of the manuscript.
Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.
Acknowledgments
Not applicable.
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Conflicts of Interest
The authors had no conflicts of interest.