A potential target gene CD63 for different degrees of intervertebral disc degeneration

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.

Low back pain (LBP) is a common cause of disability and it negatively effects on the quality of life of patients globally 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) 2 . The IDD is currently recognized among the common causes of morbidity 2,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 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 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 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 7,8 . A few studies have reported on numbers of genes and how they correlate with functional and structural changes within the IVD [9][10][11][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.  Pfirrmann classification of  IDD 13 , the gene expression profiles were obtained from 15 samples of 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  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 (GSE15227 vs. 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 14 , ECM 15 , COL2 16 , TIMP 17 , MMP 18,19 as well as ADAMTS 19,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) 21-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 1)

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. 1). Three gene modules of MEturquoise, MEblue and MEblue were obtained from the hierarchical clustering of WGCNA (Fig. 2). Among these modules, MEturquoise had 198 genes, with the most significant difference of Cor = 0.5659 and P = 0.000213 (Table 2). In the Pearson correlation analysis between DEGs and IDD, nine key genes were screened (Cor > 0.6 and P < 0.001) ( Table 3).
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. 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. 3B).
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. 5A), 3 items of molecular function (Fig. 5B) and 43 items of biological processes (Fig. 5C).
In the GSEA of CD63, 6 cellular components (|NES|> 1.8) (Table 4), 37 items of molecular function (|NES|> 1.78) ( Table 5) and 42 items of biological processes (|NES|> 1.7) ( Table 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. 6A). The molecular functions and biological processes obtained from the GSEA of CD63 were intersected with that screened using www.nature.com/scientificreports/ WGCNA in the hub modules related to IDD. Finally, one hub molecular function (Collagen binding) (Fig. 6B) and three hub biological processes (Fig. 6C) were 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. 7A). In the GSEA of CD63, 20 items of signaling pathways (Table 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. 7B) were finally obtained.

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. 8) (P < 0.01).  CD63, also known as lysosome associated membrane protein 3 (LAMP3), belongs to the transmembrane 4 superfamily (TM4SF) [24][25][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 27 .
The gene CD63 can activate the surface antigens on platelets 28 . As a sensitive marker of platelet activation, detection of that index can sense the degree of platelet activation 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 30,31 . However, there is little research on the correlation between CD63 and IDD. Through multiple gene screening modes, we finally found the hub gene   www.nature.com/scientificreports/ 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 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 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. www.nature.com/scientificreports/ 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 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 www.nature.com/scientificreports/ (Fig. 4). Previous studies had explored that TIMP and collagen were closely related to IDD, and both were important indicators of IDD 16,17 . Takawale et al 35 identified a novel mechanism in vivo for TIMP1, CD63 and www.nature.com/scientificreports/ 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  www.nature.com/scientificreports/ signaling and apoptosis inhibition 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 [36][37][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. www.nature.com/scientificreports/ www.nature.com/scientificreports/  Figure 8. 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).

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.

Data availability
The two available datasets in this paper are from the Gene Expression Omnibus (GEO) database (https:// www. ncbi. nlm. nih. gov/ geo/). www.nature.com/scientificreports/