Abstract Understanding molecular mechanisms of intervertebral disc degeneration (IDD) and providing a novel target for the treatment of IDD have important implications. We sought to explore a new promising gene target for the treatment of IDD. This study integrated 19,678 genes of 38 IDD patients from two gene datasets. Differentially Expressed Genes (DEGs) of annulus fibrosus were analyzed in groups with mild disc degeneration (MDD) and severe disc degeneration (SDD). We screened the hub gene through biological information technology (bioinformatic) methods. Then, we further validated the hub gene using annulus fibrosus and nucleus pulposus tissues from 12 patients with qRT-PCR. In addition, we explored its underlying molecular mechanism with GO, KEGG and GSEA. Through multiple screening bioinformatics methods, the hub gene CD63 was identified. The qRT-PCR explored that CD63 decreased significantly in SDD group compared to that in MDD group (P < 0.001). The GO, KEGG and GSEA of CD63 explored significant enrichment of the molecular features (P < 0.001), including the cellular component (Extracellular matrix, P < 0.001), the molecular function (collagen binding, P < 0.001), the biological processes (protein targeting, collagen fibril organization and platelet degranulation, P < 0.001) and the signaling pathways. Our research explored and validated a new regulatory gene, CD63 for different degrees of IDD. A new novel form of therapeutic target for IDD may be developed. Subject terms: Computational biology and bioinformatics, Biomarkers, Medical research Introduction Low back pain (LBP) is a common cause of disability and it negatively effects on the quality of life of patients globally^[36]1. However, the commonly reported and targeted factor for intervention is intervertebral disc degeneration (IDD). IDD plays a significant role in LBP and associates strongly with dysfunction and structural breakdown of intervertebral disc (IVD)^[37]2. The IDD is currently recognized among the common causes of morbidity^[38]2,[39]3. Nonetheless, the etiology of IDD is multifactorial and their complex mechanisms are not well understood. Some studies have reported factors including apoptosis, insufficient nutritional supply and excessive mechanical load^[40]4. Currently, treatment for IDD is largely dependent on surgical intervention with disc excision and spinal fusion, for late-stage IDD as well as symptomatic relief. In the early stages for example, conservative treatments like bed rest, painkillers, or physiotherapy are usually the preferred option^[41]5. When it seriously affects the quality of people’s life, treatment for IDD is largely dependent on surgical intervention with disc excision and spinal fusion^[42]6. It is important to note that surgical treatment does not preserve the function of disc. The pathogenesis of IDD has been vastly linked to a lot of factors for instance, spine injuries, aging, spine deformities and genes^[43]7,[44]8. A few studies have reported on numbers of genes and how they correlate with functional and structural changes within the IVD^[45]9–[46]12. Even though such studies shed light on molecular aspect of IDD, there is still little knowledge regarding gene factors and their contribution to pathogenesis of IDD as well as therapeutic targets related to the disease. Determining the onset and different degrees of intervertebral disc degeneration of IDD based on relevant regulatory genes is an important insight for the prevention of IDD and more effective treatment options for the disease. This will enable the development of novel therapies for IDD that are more specific. Therefore, it was necessary to carry out this study. It may help us further understand the molecular mechanism of IDD and provide new promising targets for the diagnosis and treatment of IDD. So, we sought to explore a new promising gene target for the treatment of IDD. Materials and METHODS Data collection and preparation We collected a total of 19,678 genes from 38 degenerative disc patients with complete clinical information from Gene Expression Omnibus (GEO) database for external analysis. Samples from two datasets [47]GSE15227 and [48]GSE23130 were selected. According to the Pfirrmann classification of IDD^[49]13, the gene expression profiles were obtained from 15 samples of [50]GSE15227 with 5 mild disc degeneration (MDD) (grade I–II) and 10 severe disc degeneration (SDD) (grade III–V). Among the 23 samples data in [51]GSE23130, gene expression profiles were obtained from 6 MDD (grade I–II) and 17 SDD (grade III-V). Both the two datasets were based on GLP1352 detection platform and come from the same microarray [U133_X3P] Affymetrix Human X3P Array. The R package “limma” was used to transform RNA expression data from a biased distribution to an approximate normal distribution. The ComBat function of package “sva” was used to remove the batch effect of samples in datasets. The gene expression profiles data was standardized with log2 logarithm and package “limma” to eliminate the heterogeneity. A total of 19,678 genes were obtained through normalization of the two datasets. Differential analysis of genes A total of 11 samples in the MDD group and 27 samples in the SDD group were analyzed. The R package “limma” was used to calculate the P-values and Fold Change (FC) of gene expression differences between the MDD group and SDD group. The P < 0.05 and |log2 FC|> 1 were selected as screening threshold for significant differential genes (DEGs). The filtered values of DEGs were extracted from data of standardized expression profiles. We also compared the similarities of the two datasets ([52]GSE15227 vs. [53]GSE23130) and performed a heatmap of the differential gene expression profiles of MDD group and SDD group in the two datasets. Weighted correlation network analysis for DEGs Weighted Correlation Network Analysis (WGCNA) of DEGs was carried out to further describe the association patterns of gene expression profiles. We calculated the Pearson correlation of DEGs with package “WGCNA”. The adjacent matrix was transformed into a topological overlap matrix. Thereafter, we calculated the dissimilarity of gene sets and gained the hierarchical clustering tree of genes with R, which were cut into different modules with a minimum number of 10 module genes. To screen out the hub module of IDD, correlations between gene modules and the grades of degenerative IVD samples were analyzed. The R package “WGCNA” was used to calculate Pearson Correlation (Cor) and P-values between DEGs and Pfirrmann classification of IDD with the screening criteria of Cor > 0.6 and P < 0.01. The key genes (Gene set A) that were closely related to IDD were identified and documented. Construction of gene co-expression network Cytoscape package “CytoNCA” were used to analyze three topological characteristics of each node in the network. This included degree of the node, number centricity and proximity to centrality. We extracted top ten common genes in each topological feature and constructed a co-expression genes (Gene set B) network with Cytoscape. Select the target gene We took an intersection of the two key gene sets that screened by WGCNA (Gene set A) and Co-expression Network (Gene set B), and having the target gene CD63. To further estimate whether CD63 was closely related to IDD, we analyzed the Pearson Correlation between CD63 and IL1^[54]14, ECM^[55]15, COL2^[56]16, TIMP^[57]17, MMP^[58]18,[59]19 as well as ADAMTS^[60]19,[61]20 which were closely related to IDD significantly. GO and KEGG analysis of gene set for hub module To further explore the mechanism of intervertebral disc degression, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG)^[62]21–[63]23 gene-set enrichment analysis of the hub modules in WGCNA. The GO defines gene function from three aspects: cellular component (CC), molecular function (MF) and biological process (BP). The KEGG provide enrichment analysis of gene pathway. The DAVID was used to perform GO and KEGG pathway enrichment analysis for the hub module. Then, we screened CC, BP, MF and signaling pathway with differential expression in the hub module of IDD. To fully understand the GO enrichment results associated with IDD, R was used to visualize the gene-set enrichment results. Gene Set enrichment analysis The GSEA were used to further identify the significantly enrichment function. The reference gene set used in this study were c5.all.v7.1.symbols.gmt and c2.cp.kegg.v7.1.symbols.gmt, with nominal P < 0.05 and FDR < 0.25 used as threshold to screen significant enrichment functions and pathways. Enrichment functions and signaling pathways of CD63 obtained in GSEA were intersected with the GO and KEGG pathways that were obtained from hub module of IDD. Patient tissue samples To further study and validate the role of CD63 in different degrees of IDD, annulus fibrosus (AF) and nucleus pulposus (NP) tissues from 12 patients (Table [64]1) were used to measure the expression of CD63 with qRT-PCR. The inclusion and exclusion criteria were as follows. Inclusion criteria: patients with mild and severe lumbar disc degeneration and treated with discectomy. Exclusion criteria: patients with intervertebral disc calcification and acute infection. The experimental protocols were approved by the Ethics Committee of Affiliated Hospital of Jining Medical University (Jining, China). Before the experiments, tissue samples of all the patients were collected with informed consent. This study performed in accordance with the declaration of Helsinki and guidelines of the Ethics Committee of Affiliated Hospital of Jining Medical University. The clinical samples were collected from patients with lumbar discectomy and were divided into MDD group (n = 6) (Supplementary Fig. [65]1) and SSD group (n = 6) (Supplementary Fig. [66]2) according to Pfirrmann classification^[67]13. In SDD of grade III-V, the boundary between annulus fibrosus and annulus fibrosus disappears, as well as the height of intervertebral space decreases, however, MDD of grade I and II only with signal intensity changes in MRI. Table 1. Clinical characteristics of patients. Item MDD group (n = 6) SDD group (n = 6) P Age (years) 33 58 55 38 39 50 40 38 41 31 35 44 Age(M ± SD) 40.50 ± 7.7395 43.17 ± 9.6833 0.6097 Clinical feature Sex (male) 4 3  > 0.999 Low back pain 5 4 Lower extremity radicular pain 6 6 [68]Open in a new tab Total RNA isolation and qRT-PCR Total RNA was extracted respectively from AF and NP tissues of 12 patients with TRIzol Reagent (Invitrogen Life Technologies, Carlsbad, CA, USA). Reverse transcribed using PrimeScript RT Master Mix (Invitrogen Life Technologies, Carlsbad, CA, USA). qRT-PCR was performed using cDNA and SYBR mixture (CWBio, Beijing, China) to quantify the mRNA expression levels of CD63 and normalized it against β-actin. Relative expression levels of mRNA were computed with the method of 2^−ΔΔ Ct. Statistical analysis The statistical analyses were calculated by GraphPad Prism 8. Significance was established at p < 0.05 and processed with GraphPad Prism 8. Ethics declarations We declare all the tissue samples of patients were collected with informed consent. The patients and their families were informed that data from the cases would be submitted for publication, and gave their consent. The experimental protocols were approved by the Ethics Committee of Affiliated Hospital of Jining Medical University (Jining, China). This study performed in accordance with the declaration of Helsinki and guidelines of the Ethics Committee of Affiliated Hospital of Jining Medical University. Results Hub module and genes Based on the threshold of P < 0.01 and |log2 FC|> 1, 290 differential genes were screened out from the two datasets. The heatmap of the expression profiles of DEGs showed the similarity of the two datasets in SDD group and MDD group (Fig. [69]1). Three gene modules of MEturquoise, MEblue and MEblue were obtained from the hierarchical clustering of WGCNA (Fig. [70]2). Among these modules, MEturquoise had 198 genes, with the most significant difference of Cor = 0.5659 and P = 0.000213 (Table [71]2). In the Pearson correlation analysis between DEGs and IDD, nine key genes were screened (Cor > 0.6 and P < 0.001) (Table [72]3). Figure 1. [73]Figure 1 [74]Open in a new tab Heatmap of differential gene expression profiles in [75]GSE15227 and [76]GSE23130. The heatmap showed the similarity of the two datasets in SDD group and MDD group. Red represents up-regulated genes and green represents down-regulated genes. Figure 2. Figure 2 [77]Open in a new tab Gene modules. Three gene modules of MEturquoise, MEblue and MEblue were obtained. The MEturquoise is the hub module. Table 2. The Pearson Correlation and P-values between gene modules and IDD. Gene module Cor P MEblue 0.4744 0.002621 MEturquoise 0.5659 0.000213 MEgrey 0.4861 0.001976 [78]Open in a new tab Table 3. The Pearson Correlation and P-values between DEGs and IDD. DEGs Cor P C2CD2 0.6501  < 0.001 GSTP1 0.6343  < 0.001 SGK1 0.6318  < 0.001 CD63 0.6222  < 0.001 LUM 0.6171  < 0.001 SCRG1 0.6131  < 0.001 C1S 0.6074  < 0.001 REEP5 0.6034  < 0.001 CLIC1 0.6016  < 0.001 [79]Open in a new tab Gene co-expression network Through the analysis of WGCNA, a total of 81 gene nodes and 43 interaction pairs were obtained with the threshold of Cor > 0.6 as explored in the Gene Co-expression Network (Fig. [80]3A). It was found that the top 10 genes related to IDD in each Topological Features (Degrees, Meso-centricity, Proximity-centricity) including CD63, PAM, SSR4 and RPS19 (Fig. [81]3B). Figure 3. [82]Figure 3 [83]Open in a new tab Gene co-expression network and topological features. (A) A total of 81 gene nodes and 43 interaction pairs obtained in the Weighted Correlation Network Analysis with the threshold of Cor > 0.6. (B) Intersection of the top 10 genes in each Topological Features (Betweenness, Closeness and Degree) related to intervertebral disc degeneration. Selection of target genes The genes that closely related to IDD, based on WGCNA (Gene set A) and Co-expression Network (Gene set B) screening were intersected and the hub gene CD63 obtained (Fig. [84]4A). Pearson Correlation Analysis between CD63 and IL1, ECM, COL2A1, TIMP, MMP and ADAMTS were performed (Fig. [85]4B). The results explored that CD63 was significantly related to TIMP1 (Cor = 0.789), TIMP2 (Cor = 0.887), TIMP3 (Cor = 0.856) and COL2A1 (Cor = 0.849) (Fig. [86]4C). Figure 4. [87]Figure 4 [88]Open in a new tab Selection of target genes. (A) The genes that closely related to intervertebral disc degeneration screened in weighted correlation network analysis (Gene set A) and co-expression network (Gene set B). (B) Pearson correlation analysis between CD63 and IL1, ECM, COL2A1, TIMP, MMP and ADAMTS. (C) Significant relatedness of CD63 with TIMP1, TIMP2, TIMP3 and COL2A1. GO, KEGG and GSEA The analysis of GO (CC, MF and BP) enrichment for IDD was carried out. Genes in MEturquoise, the hub module associated with IDD, was significantly enriched in 7 cellular components (Fig. [89]5A), 3 items of molecular function (Fig. [90]5B) and 43 items of biological processes (Fig. [91]5C). Figure 5. [92]Figure 5 [93]Open in a new tab The analysis of GO (Cellular Component, Molecular Function and Biological Process) enrichment for intervertebral disc degeneration. (A) Genes in the hub module associated with intervertebral disc degeneration significantly enriched in 7 cellular components, (B) 3 items of molecular function and (C) 43 items of biological processes. In the GSEA of CD63, 6 cellular components (|NES|> 1.8) (Table [94]4), 37 items of molecular function (|NES|> 1.78) (Table [95]5) and 42 items of biological processes (|NES|> 1.7) (Table [96]6) were obtained. Cellular components related to IDD that was screened from GSEA were intersected with that screened using WGCNA in hub module. The hub cellular component (Extracellular matrix) was finally obtained (Fig. [97]6A). The molecular functions and biological processes obtained from the GSEA of CD63 were intersected with that screened using WGCNA in the hub modules related to IDD. Finally, one hub molecular function (Collagen binding) (Fig. [98]6B) and three hub biological processes (Fig. [99]6C) were obtained. Table 4. The GSEA of CD63 for cellular components. Cellular components P FDR NES Extracellular exosome  < 0.001 0.0206 1.9779 Extracellular matrix  < 0.001 0.0932 1.8534 Cytosol  < 0.001 0.0777 1.8259 Extracellular space 0.0020 0.0641 1.8163 Membrane  < 0.001 0.0528 1.8147 Focal adhesion 0.0020 0.0529 1.8016 [100]Open in a new tab Table 5. The GSEA of CD63 for molecular function. Term P FDR NES Extracellular matrix structural constituent  < 0.001 0.0654 1.9780 Translation elongation factor activity  < 0.001 0.0327 1.9761 WW domain binding  < 0.001 0.0232 1.9741 Fibroblast growth factor binding  < 0.001 0.0225 1.9589 Growth factor binding  < 0.001 0.0195 1.9501 Glycosaminoglycan binding  < 0.001 0.0383 1.8907 Protease binding  < 0.001 0.0397 1.8764 Disordered domain specific binding  < 0.001 0.0359 1.8720 Calcium dependent phospholipid binding 0.0021 0.0323 1.8695 RNA helicase activity  < 0.001 0.0372 1.8513 Beta catenin binding  < 0.001 0.0372 1.8452 Translation regulator activity  < 0.001 0.0359 1.8410 Cell adhesion molecule binding  < 0.001 0.0369 1.8351 SMAD binding  < 0.001 0.0347 1.8341 GTPase activity  < 0.001 0.0324 1.8335 Extracellular matrix structural constituent conferring compression resistance 0.0020 0.0306 1.8327 Ubiquitin like protein ligase binding  < 0.001 0.0297 1.8314 Kinase regulator activity  < 0.001 0.0293 1.8289 Copper ion binding  < 0.001 0.0312 1.8229 Translation regulator activity nucleic acid binding  < 0.001 0.0298 1.8217 Calcium dependent protein binding  < 0.001 0.0296 1.8196 ATPase binding  < 0.001 0.0283 1.8191 GDP binding  < 0.001 0.0262 1.8181 Ribosome binding  < 0.001 0.0274 1.8181 Hyaluronic acid binding 0.0021 0.0255 1.8176 Cadherin binding  < 0.001 0.0258 1.8137 Unfolded protein binding  < 0.001 0.0263 1.8082 Monosaccharide binding  < 0.001 0.0258 1.8080 mRNA binding  < 0.001 0.0252 1.8070 Heparin binding  < 0.001 0.0244 1.8066 Nucleobase containing compound transmembrane transporter activity  < 0.001 0.0252 1.8013 mRNA 3 UTR binding  < 0.001 0.0248 1.8002 Metalloendopeptidase inhibitor activity 0.0082 0.0255 1.7952 L ascorbic acid binding  < 0.001 0.0271 1.7904 Heat shock protein binding  < 0.001 0.0274 1.7877 Extracellular matrix structural constituent conferring tensile strength 0.0021 0.0270 1.7864 Collagen binding 0.0019 0.0274 1.7813 [101]Open in a new tab Table 6. The GSEA of CD63 for biological processes. Term P FDR NES Keratan sulfate metabolic process  < 0.001 0.2063 1.9954 Cartilage development 0.0020 0.2377 1.7909 Protein targeting to membrane 0.0021 0.2327 1.7866 Pathway restricted SMAD protein phosphorylation 0.0020 0.2220 1.7857 Entrainment of circadian clock  < 0.001 0.2138 1.7840 Relaxation of cardiac muscle  < 0.001 0.2035 1.7837 Cellular response to gamma radiation  < 0.001 0.1974 1.7820 Chondroitin sulfate biosynthetic process  < 0.001 0.1929 1.7799 Proteoglycan biosynthetic process 0.0020 0.1965 1.7742 Artery development  < 0.001 0.1943 1.7705 Connective tissue development  < 0.001 0.1991 1.7647 Chondroitin sulfate proteoglycan biosynthetic process 0.0020 0.1938 1.7636 Negative regulation of cellular response to growth factor stimulus  < 0.001 0.1867 1.7635 RNA polyadenylation  < 0.001 0.1813 1.7630 Bone growth  < 0.001 0.1820 1.7599 Chondrocyte differentiation 0.0020 0.1785 1.7586 Aminoglycan metabolic process 0.0020 0.1850 1.7522 Labyrinthine layer blood vessel development 0.0021 0.1821 1.7504 Protein targeting  < 0.001 0.1774 1.7501 Regulation of morphogenesis of a branching structure  < 0.001 0.1742 1.7492 Positive regulation of morphogenesis of an epithelium  < 0.001 0.1709 1.7487 Growth plate cartilage morphogenesis  < 0.001 0.1794 1.7421 Negative regulation of transmembrane receptor protein serine threonine kinase signaling pathway  < 0.001 0.1748 1.7419 Sulfur compound biosynthetic process  < 0.001 0.1711 1.7413 Entrainment of circadian clock by photoperiod  < 0.001 0.1675 1.7409 Chondrocyte development 0.0020 0.1633 1.7408 Extracellular structure organization  < 0.001 0.1637 1.7378 Ovulation  < 0.001 0.1618 1.7369 Chaperone mediated protein folding  < 0.001 0.1667 1.7316 Regulation of protein maturation 0.0041 0.1630 1.7316 Sulfur compound catabolic process 0.0020 0.1678 1.7279 Chondroitin sulfate proteoglycan metabolic process 0.0020 0.1679 1.7260 Cell recognition  < 0.001 0.1694 1.7231 Bone development  < 0.001 0.1663 1.7229 Collagen fibril organization 0.0042 0.1640 1.7223 Regulation of DNA templated transcription in response to stress  < 0.001 0.1609 1.7222 Cellular aldehyde metabolic process 0.0021 0.1628 1.7191 Mesenchyme morphogenesis  < 0.001 0.1680 1.7146 Platelet degranulation 0.0020 0.1718 1.7114 Cranial nerve development  < 0.001 0.1783 1.7064 Aminoglycan catabolic process 0.0063 0.1782 1.7044 Regulation of receptor biosynthetic process 0.0021 0.1783 1.7025 [102]Open in a new tab Figure 6. [103]Figure 6 [104]Open in a new tab GO (cellular component, molecular function and biological process) related to intervertebral disc degeneration in the hub modules intersected with gene set enrichment analysis of CD63. The hub cellular component (extracellular matrix) (A) molecular function (collagen binding) (B) as well as three biological processes (protein targeting, collagen fibril organization and Platelet degranulation) (C) were finally obtained. Signaling pathways The enrichment analysis of KEGG pathways were performed for the genes in hub module (MEturquoise) that were screened using WGCNA. Results showed that the genes were significantly enriched in 12 signaling pathways (Fig. [105]7A). In the GSEA of CD63, 20 items of signaling pathways (Table [106]7) (|NES|> 1.55) were enriched. The signaling pathways obtained from GSEA of CD63 were intersected with genes screened using WGCNA in the hub module related to IDD, and 4 hub signaling pathways (Fig. [107]7B) were finally obtained. Figure 7. [108]Figure 7 [109]Open in a new tab KEGG pathways related to intervertebral disc degeneration in the hub modules intersected with gene set enrichment analysis of CD63. (A) KEGG pathways significantly enriched in 12 signaling pathways. (B) 4 hub pathways were finally obtained when intersected KEGG pathways with gene set enrichment analysis of CD63. Table 7. The GSEA of CD63 for signaling pathways. Term P FDR NES Cardiac muscle contraction 0.0079 0.0641 1.8305 Alzheimers disease  < 0.001 0.0676 1.7623 Huntingtons disease  < 0.001 0.0576 1.7500 Glycolysis gluconeogenesis 0.0039 0.0798 1.7115 Oxidative phosphorylation  < 0.001 0.0673 1.7080 Parkinsons disease 0.0020 0.0809 1.6816 Galactose metabolism 0.0078 0.0722 1.6767 Glycosaminoglycan biosynthesis chondroitin sulfate 0.0020 0.0801 1.6585 Selenoamino acid metabolism 0.0040 0.0834 1.6444 ECM receptor interaction 0.0020 0.0849 1.6342 Arrhythmogenic right ventricular cardiomyopathy  < 0.001 0.0873 1.6223 Spliceosome  < 0.001 0.0912 1.6101 Ribosome 0.0061 0.1003 1.5917 Glutathione metabolism 0.0059 0.1057 1.5809 N glycan biosynthesis 0.0040 0.1006 1.5790 Purine metabolism 0.0060 0.0948 1.5783 Focal adhesion 0.0059 0.0974 1.5704 Proteasome  < 0.001 0.1037 1.5548 Gap junction 0.0020 0.1033 1.5501 Notch signaling pathway 0.0079 0.0984 1.5501 [110]Open in a new tab Expression of CD63 in human AF and NP tissues To further study and validate the role of CD63 in different degrees of IDD, we measured the expression of CD63 in AF and NP tissues from MDD patients (n = 6) and SDD patients (n = 6) with qRT-PCR. Results showed that the mRNA expression of CD63 in AF and NP tissues was markedly downregulated in SDD group compared to that in MDD group (Fig. [111]8) (P < 0.01). Figure 8. [112]Figure 8 [113]Open in a new tab The expression of CD63 in annulus fibrosus and nucleus pulposus tissues from MDD (n = 6) and SDD (n = 6) patients with qRT-PCR. The mRNA expression of CD63 in annulus fibrosus and nucleus pulposus was markedly downregulated in SDD group compared to that in MDD group (P < 0.01). Discussion This study integrated 19,678 genes of 38 IDD patients from two gene datasets. We carried on multiple gene screening modes (DEGs Analysis, Pearson Correlation Analysis, WGCNA and Topological Characteristics) via biological information technology approach. Finally, the new regulatory gene for IDD, CD63, was identified. CD63, also known as lysosome associated membrane protein 3 (LAMP3), belongs to the transmembrane 4 superfamily (TM4SF)^[114]24–[115]26. The TM4SF members are related to each other and form a huge TM4SF network with some extra family proteins, which play vital roles in molecular metabolism^[116]27. The gene CD63 can activate the surface antigens on platelets^[117]28. As a sensitive marker of platelet activation, detection of that index can sense the degree of platelet activation^[118]29, and then make a timely diagnosis and curative measure against the disease. Previous studies have reported the positive role of CD63 in the suppression of melanoma whereby it acted as a vital sign in patient assessment^[119]30,[120]31. However, there is little research on the correlation between CD63 and IDD. Through multiple gene screening modes, we finally found the hub gene CD63. In the datasets, it was a significantly downregulated gene in SDD group compared to that in MDD group. Finally, the results of qRT-PCR of annulus fibrosus and nucleus pulposus tissues further validated that the expression of CD63 was markedly downregulated in SDD group. Different degrees of IDD may be closely related to the reduced expression of CD63. It may become a new promising gene for IDD and help us further enrich the therapeutic targets of IDD. Extracellular matrix (ECM) located outside one or more cells and provided structural support. The main components of ECM are proteoglycans and collagen II that maintain the physiological function and stability of IVD. The components are reduced during the process of IDD, thus can be used as an important sign for IDD^[121]32. The ECM plays critical role in the maintenance of steady state for IVD and different degrees of IDD. CD63 was a cell surface binding partner for Tissue Inhibitor of Metalloproteinases-1 (TIMP-1) which play important role in the degradation and synthesis of matrix^[122]33. The result of GO enrichment analysis also showed that CD63 play an important role in cell matrix. This is consistent with previous studies. It suggested that CD63 may participate in the regulation of IDD through ECM. The screened MF of CD63 obtained using the conjoint analysis of WGCNA and GSEA were "collagen binding". Collagen, a main component of IVD, was fully demonstrated in IDD^[123]34. Collagen binding implies that a group of fibrous proteins with high tensile strength interact with collagen selectively and non-covalently. The Pearson Correlation Analysis further indicated that CD63 was significantly related with TIMP and collagen (Fig. [124]4). Previous studies had explored that TIMP and collagen were closely related to IDD, and both were important indicators of IDD^[125]16,[126]17. Takawale et al^[127]35 identified a novel mechanism in vivo for TIMP1, CD63 and collagen synthesis, as well as CD63 shown to exist as a cell surface receptor for TIMP1. Their study also showed that TIMP1 mediates an association between CD63 and integrin β1, leading to de novo collagen synthesis on cardiac fibroblasts. The screened MF of CD63 (Collagen binding) explored potential molecular mechanism of CD63, TIMP1 and collagen in different degrees of IDD. Using yeast two-hybrid screening, Jung et al. identified CD63 as a cell surface binding partner for TIMP1 which played a critical role in TIMP1-mediated cell survival signaling and apoptosis inhibition^[128]33. In the future, the relationship between CD63 and TIMP1 would be reported in more diseases. In addition, the significant enrichment of “protein targeting” and “collagen fibril organization” in BP, as well as “proteasome” and “oxidative phosphorylation” in the pathway of CD63 explored new basis and implications for the study of CD63 in IDD. Many previous studies also have reported the causes of IDD for example, decreased expression of collagen II and proteoglycan, increased expression of enzymic degradation in ECM and apoptosis in nucleus pulposus cells^[129]36–[130]38. Although we have tested the predictive results through clinical samples and explored their possible mechanism, the sample sizes may be still insufficient according to the difference between mild and severe disc degeneration groups studied with bioinformatic analysis. To the best of our knowledge, similar reports are rare. Therefore, in prospective studies, we intend to analyze the underlying molecular mechanisms mechanism of CD63 for IDD with multiple biological methods. Conclusion Our study explored and validated a new target gene CD63 for different degrees of IDD. The study and its findings are important as they form a basis upon which new therapeutic targets for IDD can be identified. Supplementary Information [131]Supplementary Information.^ (690.8KB, docx) Acknowledgements