Abstract
Aims
This study aimed, through bioinformatics analysis and in vitro experiment validation, to identify the key extracellular proteins of intervertebral disc degeneration (IDD).
Methods
The gene expression profile of GSE23130 was downloaded from the Gene Expression Omnibus (GEO) database. Extracellular protein-differentially expressed genes (EP-DEGs) were screened by protein annotation databases, and we used Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) to analyze the functions and pathways of EP-DEGs. STRING and Cytoscape were used to construct protein-protein interaction (PPI) networks and identify hub EP-DEGs. NetworkAnalyst was used to analyze transcription factors (TFs) and microRNAs (miRNAs) that regulate hub EP-DEGs. A search of the Drug Signatures Database (DSigDB) for hub EP-DEGs revealed multiple drug molecules and drug-target interactions.
Results
A total of 56 EP-DEGs were identified in the differential expression analysis. EP-DEGs were enriched in the extracellular structure organization, ageing, collagen-activated signalling pathway, PI3K-Akt signalling pathway, and AGE-RAGE signalling pathway. PPI network analysis showed that the top ten hub EP-DEGs are closely related to IDD. Correlation analysis also demonstrated a significant correlation between the ten hub EP-DEGs (p＜0.05), which were selected to construct TF–gene interaction and TF–miRNA coregulatory networks. In addition, ten candidate drugs were screened for the treatment of IDD.
Conclusion
The findings clarify the roles of extracellular proteins in IDD and highlight their potential as promising novel therapeutic targets.
Cite this article: Bone Joint Res 2023;12(9):522–535.
Article focus
-
What are the extracellular protein expression changes in intervertebral disc degeneration (IDD) processes?
-
What are the function and pathways of extracellular proteins in IDD processes?
-
What is the relationship between transcription factors (TFs), microRNAs (miRNAs), and candidate drugs regulating extracellular proteins in IDD processes?
Key messages
-
Extracellular proteins mediate the progression of IDD.
-
TFs, miRNAs, and candidate drugs regulating extracellular proteins are of great significance in the development of IDD.
Strengths and limitations
-
The present study revealed that extracellular proteins are a potential therapeutic target for IDD.
-
Experimental verification and a large-scale sample size are still needed for better reliability of the results in the future.
Introduction
Intervertebral disc degeneration (IDD) is a common disease in orthopaedics and one of the leading causes of chronic lower back pain (LBP), which seriously reduces patients’ quality of life and represents a sizeable economic burden.1-3 It is estimated that approximately 35% of subjects aged between 20 and 39 years experience at least one IDD at the lumbar spine level, while subjects aged between 60 and 80 years have IDD of varying degrees.4 The current treatment for IDD is mainly pain relief. For early treatment, conservative methods such as drugs and physical therapy are included, while surgical treatment is the main option for advanced IDD patients, including disc arthroplasty and spinal fusion.5,6 In recent decades, with the rapid development of molecular biology techniques, research on the molecular mechanism of IDD has progressed, but no effective treatment to delay or reverse the chronic progression of IDD exists. Therefore, understanding the molecular mechanism of IDD and exploring the potential therapeutic targets are essential.
A healthy intervertebral disc (IVD) consists of three parts, the central nucleus pulposus, the peripheral annulus fibrosus, and the cartilage endplates on both sides, which provide the structural support for movement and spinal load.7,8 Many factors, such as ageing, mechanical stress, smoking, obesity, trauma, and genetics, lead to reduced IVD height, fissure or rupture of the annulus fibrosus, loss of nucleus pulposus collagen and water, and calcification of the cartilage endplates. Subsequently, pathological changes such as bulging, protrusion, and prolapse occur, and the clinical manifestations of LBP appear in the corresponding parts.9-11 In recent years, increasing evidence has suggested that some secreted proteins in the extracellular environment of IVD, namely extracellular matrix (ECM) proteins, play an important role in maintaining the structure and function of IVD. Therefore, studies of extracellular proteins in IDD processes can lead to the discovery of new therapeutic targets to delay or reverse IDD.12-14
Currently, nucleus pulposus cells show the highest degree of remodelling and maintain ECM homeostasis during IDD. Under normal conditions, ECM anabolism and catabolism are in homeostasis. IDD typically occurs when the homeostatic balance of the ECM is disturbed by various stimuli.8,12 With the chronic progression of IDD, proteoglycan and type II collagen content notably decreases, whereas type I collagen content notably increases, resulting in reduced water absorption, a dysfunctional mechanical microenvironment, and reduced resistance to ECM loading.15 Furthermore, the aberrant expression of various proteases in degenerated disc tissue and their roles in ECM proteolysis have been reported, representing potential biomarkers or therapeutic targets for IDD.12,16
Bioinformatics analysis of microarray transcription profiles is a novel approach for exploring disease pathogenesis, identifying disease biomarkers, and discovering therapeutic targets.4,16 The Gene Expression Omnibus (GEO) database contains microarray expression profiling data for many diseases.17 We downloaded the IDD patient IVD tissue microarray gene expression profiling dataset (GSE23130) from the GEO database and used R software to screen for differentially expressed genes (DEGs) between IVD tissue samples of early and late IDD. Next, we examined extracellular protein-DEGs (EP-DEGs). Subsequently, we performed functional enrichment and pathway enrichment analysis of EP-DEGs using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. A protein-protein interaction (PPI) network of EP-DEGs was also established to screen functional modules and hub genes. Finally, we collected nucleus pulposus samples from IDD patients for histological evaluation and validated the expression of hub EP-DEGs using quantitative real-time polymerase chain reaction (qRT-PCR) analysis (Figure 1). Our findings clarify the roles of extracellular proteins in IDD and highlight their potential as promising novel therapeutic targets.
Fig. 1
Flowchart of the identification of the key extracellular proteins as potential targets in intervertebral disc degeneration. DEG, differentially expressed gene; EP-DEG, extracellular protein-DEG; GO, Gene Ontology; H&E, haematoxylin and eosin; HPA, Human Protein Atlas; IDD, intervertebral disc degeneration; miRNA, microRNA; PPI, protein-protein interaction; qRT-PCR, quantitative real-time polymerase chain reaction; TF, transcription factor.
Methods
Data acquisition and processing
The GEO database is an open data access platform.18 It contains the most extensive collection of microarray and high-throughput sequencing gene expression profiling data to date. We used R software (ver. 3.6.3; R Foundation for Statistical Computing, Austria) and the GEOquery package to download and analyze the reliable GSE23130 IDD expression profile dataset from the GEO database.19 The dataset source is human, and the platform is based on the GPL1352 [U133_X3P] Affymetrix Human X3P Array (Affymetrix, USA). The GSE23130 dataset includes 23 IVD tissue samples with different degrees of degeneration. According to the Pfirrmann grading system,20 the dataset contained one case of grade I, five cases of grade II, nine cases of grade III, five cases of grade IV, and three cases of grade V. We defined grades I to III as early IDD and grades IV and V as late IDD in the subsequent analyses. We applied the affy package to the original CEL file of the GSE25101 dataset through the robust multiarray average (RMA) algorithm for background correction, data normalization, and inter-sample correction, after which we applied the corresponding GPL1352 platform in Bioconductor.21 The gene annotation file annotates the probe matrix.
Screening of DEGs and extracellular protein-DEGs
After data pre-processing, DEGs were screened out by the limma package22 using |log2 fold change (FC)| ≥ 1, p < 0.05 (independent-samples t-test) as the differential cut-off standard. Next, we used the Human Protein Atlas (HPA) protein annotation database to download the extracellular protein gene list,23 and the UniProt database to download the extracellular protein gene list GO:0005576.24-26 Finally, we combined the results with DEGs to screen out EP-DEGs and then analyzed the differential expression of EP-DEGs in the early and late IDD groups.
Functional enrichment and pathway analysis of EP-DEGs
We performed GO and KEGG pathway enrichment analysis on EP-DEGs using the ClusterProfiler package27 in Bioconductor.21 We selected p < 0.05 (independent-samples t-test) and counts ≥ 2 as cut-off values, and finally obtained the biological process (BP), cellular component (CC), and molecular function (MF) of the EP-DEGs. Before performing KEGG pathway enrichment analysis, we divided the EP-DEGs into up- and down-regulated groups. We identified pathways enriched by EP-DEGs using p < 0.05 (independent-samples t-test) and counts ≥ 2 as cut-off values. Finally, we presented the enrichment analysis results using the GOplot and ggplot2 packages in R software (ver. 3.2.3; R Foundation for Statistical Computing).
Protein-protein interaction network construction and identification of hub genes
The STRING database was used to evaluate PPIs in functional protein association networks.28,29 We constructed the PPI network of EP-DEGs using STRING (v. 3.9.0) and set the biological species as “Homo sapiens”. The reliability was set at > 0.4, and then we obtained the PPI network data and imported the information obtained into the Cytoscape software for interaction network mapping and further analysis.30 The molecular complex detection (MCODE) was used to identify the key modules from the PPI network. The MCC method is the most accurate in cytoHubba.31 We used the cytoHubba plugin in Cytoscape to screen out the top ten genes as hub genes. Finally, we used the ggplot2 package in R software (ver. 3.6.3) for the statistical analysis and visualization of the hub gene expression and correlation of both groups.
Transcription factor–gene interactions
NetworkAnalyst,32 an online visual analysis platform for gene expression analysis and meta-analysis, was used to identify TF–gene interactions with the identified genes. The TF–gene interaction network was generated by the JASPAR database included in the NetworkAnalyst platform.33
TF–miRNA coregulatory network
TF–microRNA (miRNA) coregulatory interactions were collected from the RegNetwork repository, which facilitates the detection of miRNAs and regulatory TFs that regulate DEGs of interest at post-transcriptional and transcriptional levels. NetworkAnalyst is used to navigate datasets in order to identify biological features and functions that lead to valid biological hypotheses; in this study, NetworkAnalyst was used to visualize the TF–miRNA coregulatory network.
Identification of candidate drugs
Identification of protein–drug interactions is an important aspect in the screening of potential drugs for IDD. Based on our screened hub EP-DEGs, drug molecules were identified using the Drug Signatures Database (DSigDB), which was accessed through Enrichr,34 a popular analysis platform with many gene set libraries for exploring genome-wide gene set enrichment.
Clinical nucleus pulposus surgical sample collection
MRI was conducted as part of patient care. According to the Pfirrmann grading system standard,15,20 six cases of grade II or III nucleus pulposus tissue specimens were defined as the early IDD group (n = 6; two males; four females; aged 16 to 68 years (mean 38.17 years (standard deviation (SD) 18.54)), and six cases of grade IV or V nucleus pulposus tissue specimens were defined as the late IDD group (n = 6; three males; three females; aged 44 to 60 years (mean 53.33 years (SD 5.96)). None of the patients included in the study had a history of tumour, tuberculosis, or other infectious diseases. All the clinical tissue samples used in this study were obtained with the patients’ consent and signed informed consent was obtained. The Ethics Committee of Lanzhou University Second Hospital approved all the procedures.
Haematoxylin and eosin staining
The human nucleus pulposus tissue specimens were fixed in 4% paraformaldehyde and were then routinely embedded in paraffin before being cut into 5 μm thick slices. Haematoxylin and eosin (H&E) staining was performed using a H&E Staining Kit (Solarbio, China) according to the manufacturer’s instructions. Briefly, all sections were deparaffinized, rehydrated, and stained in haematoxylin solution for five to 20 minutes. After differentiation in 1% acidic alcohol for 30 seconds, the sections were stained in eosin solution for one to two minutes and washed. Finally, the slides were mounted with neutral balsam and observed under an optical microscope (Olympus, Japan) for image processing.
Safranin O–Fast Green staining
The nucleus pulposus tissue samples were fixed with paraformaldehyde, dehydrated, and then embedded in paraffin. Staining was performed using the Safranin O–Fast Green Cartilage Staining Kit (Solarbio) according to the manufacturer’s instructions. Briefly, after deparaffinization of the nucleus pulposus tissue samples, freshly prepared Weigert dye solution was added for three to five minutes, followed by washing with water. The slides were differentiated in acidic differentiation solution for 15 seconds, washed with distilled water for ten minutes, soaked in Fast Green staining solution for five minutes, and then briefly washed with a weak acid solution for ten to 15 seconds to remove any excess Fast Green. The samples were then introduced into the Safranin O stain for five minutes, dehydrated using 95% ethanol and absolute ethanol, and cleared using xylene. Finally, the slides were mounted with neutral balsam and observed under an optical microscope (Olympus) for image processing.
Toluidine blue staining
The nucleus pulposus tissue samples were fixed with paraformaldehyde, dehydrated, and then embedded in paraffin. Staining was performed according to the instructions for the Toluidine blue staining solution (the 1% phosphate method). Briefly, the nucleus pulposus tissue samples were deparaffinized, stained with 1% Toluidine blue staining solution for two hours, and washed with water. Finally, the slides were mounted with neutral balsam and observed under an optical microscope (Olympus) for image processing.
Total RNA isolation and quantitative real-time polymerase chain reaction analysis
We extracted total RNA from the nucleus pulposus tissues of 12 patients separately with TRIzol reagent (Invitrogen; Thermo Fisher Scientific, USA). Reverse transcription was performed using PrimeScript RT Master Mix (Takara Bio, Japan). RNA concentration was measured on a NanoDrop to determine the absorbance at 260 nm and 280 nm. An A260/A280 ratio between 1.8 and 2.0 suggests that the RNA quality is up to standard and can be used for subsequent experiments. According to the manufacturer’s instructions, reverse transcription was performed using a cDNA reverse transcription kit (Takara, China), which synthesized complementary DNA (cDNA) from 1 μg of total RNA. The messenger RNA (mRNA) level was determined by qRT-PCR according to standard procedures using the Light Cycler 480 System (Roche, Germany). We used a 25 μl reaction system, including 12.5 μl of TB Green Premix Ex Taq II, 1 μl of polymerase chain reaction (PCR) forward primer, 1 μl of PCR reverse primer, 2 μl of cDNA solution, and 8.5 μl of sterile water. The PCR reaction conditions were as follows: pre-denaturation (95°C for 30 seconds, one cycle), PCR reaction (95°C for 5 seconds, 60°C for 30 seconds, 40 cycles), dissolution (95°C for 5 seconds, 60°C for one minute, one cycle), and cool-down (50°C for 30 seconds, one cycle). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used for normalization and the results detected by qRT-PCR were recorded as Ct values and analyzed using the 2-ΔΔCt method. The primer sequences are listed in Table I.
Table I.
Primer sequences used for quantitative real-time polymerase chain reaction (species: Homo sapiens).
| Gene name | Primer sequences (5′-3′) |
|---|---|
| COL1A1 | F: 5'-CACACGTCTCGGTCATGGTA -3' R: 5'-AAGAGGAAGGCCAAGTCGAG-3' |
| COL1A2 | F: 5'- GAGGGCAACAGCAGGTTCACTTA-3' R: 5'- TCAGCACCACCGATGTCCAA-3' |
| COL3A1 | F: 5'-CCAGGAGCTAACGGTCTCAG-3' R: 5'-CAGGGTTTCCATCTCTTCCA-3' |
| COL4A2 | F: 5′-TGGGACAGACGAGACAACAG -3′ R: 5′-CAACGGTATTTGGGAGAACAT -3′ |
| COL4A1 | F: 5'-CTAATGTCACAACATGGTGCTAC-3' R: 5'-GCAGGGTGTGTTAGTTACGC-3' |
| COL6A3 | F: 5′-ATGAGGAAACATCGGCACTTG-3′ R: 5′-GGGCATGAGTTGTAGGAAAGC-3′ |
| COL2A1 | F: 5'- GGGATCGTGGTGACAAAGGT-3' R: 5'- CTGGGCAGCAAAGTTTCCAC -3' |
| LUM | F: 5'-CCTGGAGGTCAATCAACTTGAGA-3' R: 5'-CGATTGCCATCCAAACGCAA-3' |
| FN1 | F: 5'-TGGAGATGAGTGGGAACGAA -3' R: 5'-CGGTCCCACTTCTCTCCAAT -3' |
| COL15A1 | F: 5′-CAACGGCTGCAATTGGTCTC-3′ R: 5′-AAGACCTTGATGTGCCCAGG-3′ |
| GAPDH | F: 5′-CCGCATCTTCTTTTGCGTC-3′ R: 5′-CCCAATACGACCAAATCCGTTG-3′ |
-
COL1A1, collagen type I alpha 1 chain; F, forward; FN1, fibronectin 1; GAPDH, glyceraldehyde 3-phosphate dehydrogenase; LUM, lumican; R, reverse.
Statistical analysis
Public gene expression data were analyzed using the GEOquery package of R software, the limma package, the ClusterProfiler package, and Cytoscape software (ver. 3.9.0). Statistical analysis of experimental data was performed using GraphPad Prism (ver. 8.0; GraphPad Software, USA). Measurement data are described as means and SDs (normal distribution). We used independent-samples t-test or Wilcoxon signed-rank test for the statistical analysis of two paired groups. Pearson’s correlation analysis was used to explore correlations between the expression levels of two paired groups. Statistical significance was set at p < 0.05.
Results
Identification of DEGs
We performed differential gene expression analysis of the early and late IDD groups to compare the differences in gene expression. First, the GSE23130 dataset was standardized, and the results after data correction were displayed as box plots (Figure 2a). The mean gene expression of each sample in the dataset was almost identical, indicating that the data were successfully normalized, and the consistency was high. A total of 428 DEGs were screened, among which 400 were highly expressed DEGs and 28 were low-expressed DEGs (Figure 2b). Subsequently, the top 20 DEGs with high and low expression are displayed in heat maps, and these genes differed significantly between the early and late IDD groups. The top three upregulated genes with the most apparent differences were lumican (LUM), fibronectin 1 (FN1), and collagen type III alpha 1 chain (COL3A1). The top three downregulated genes were RAB11B, GSK3A, and CNOT3 (Figure 2c).
Fig. 2
Identification of differentially expressed genes (DEGs). a) GSM numbers correspond to Gene Expression Omnibus (GEO) samples. Box plots of gene probe expression levels in the samples. The mean gene expression values in both samples are almost identical, indicating normalized data with high consistency. b) Volcano plot of DEGs in the early and late intervertebral disc degeneration (IDD) groups. Red dots indicate upregulation, and blue dots indicate downregulation. c) Heat map of the top 20 DEGs for high and low expression. Red represents upregulated genes, and blue represents downregulated genes. *Independent-samples t-test.
Screening of extracellular protein-DEGs
We screened EP-DEGs from DEGs and analyzed them with reference to the extracellular protein genes annotated in existing public databases to screen out the DEGs that encode the extracellular proteins in the early and late IDD groups. The genes encoding extracellular proteins marked in the HPA and Uniprot databases were intersected with DEGs. A total of 56 EP-DEGs were screened out, including 54 upregulated and two downregulated genes (Figure 3a, Supplementary Table i). The top ten upregulated genes with the most apparent differences in the late IDD group were LUM, FN1, COL3A1, COL1A2, HTRA1, MGP, ASPN, COL1A1, MXRA5, and B2M, and the two downregulated genes were GSN and DMKN (Figure 3b). Expression changes and cluster analysis of the 56 EP-DEGs are shown as a heat map (Figure 3c).
Fig. 3
Screening of extracellular protein-differentially expressed genes (EP-DEGs). a) Venn diagram showing 56 EP-DEGs screened by intersecting the genes encoding extracellular proteins and DEGs downloaded from the Human Protein Atlas (HPA) and Uniprot databases. b) Volcano plot of EP-DEGs in the early and late intervertebral disc degeneration (IDD) groups. The red dots indicate upregulation, and the blue dots indicate downregulation. c) Heat map of 56 EP-DEGs with high and low expression. Red represents upregulated, and blue represents downregulated genes. COL1A1, collagen type I alpha 1 chain; LUM, lumican. *Independent-samples t-test.
Functional enrichment and pathway analysis of EP-DEGs
We performed GO and KEGG enrichment analyses on EP-DEGs to investigate their biological functions. The GO terms are composed of three parts: CC, BP, and MF. The BP DEGs were involved in extracellular structure organization, ECM organization, ageing, collagen-activated signalling pathway, and bone growth. CC analysis revealed that EP-DEGs were markedly enriched in the collagen-containing ECM, endoplasmic reticulum lumen, basement membrane, Golgi lumen, and lysosomal lumen. MF analysis revealed significantly enriched terms for ECM structural constituent, collagen binding, growth factor binding, and platelet-derived growth factor binding (Supplementary Figure aa). We used the cnetplot function in the ClusterProfiler package to display the enriched EP-DEGs in BP, CC, and MF (Figures 4a to 4c, Table II). Among them, AEBP1, ANX2, COL1A1, COL1A2, COL2A1, COL3A1, COL4A1, COL4A2, COL6A3, COL8A2, COL15A1, FN1, and LUM are enriched in at least two aspects of BPs, CCs, and MFs (Figures 4a to 4c).
Fig. 4
Circle graph of the functional enrichment analysis of extracellular protein-differentially expressed genes (EP-DEGs). a) to c) The circle graph shows the EP-DEGs enriched in the Gene Ontology (GO) categories (BP, CC, and MF). The blue points represent the GO categories; the size of a point indicates the number of genes included. d) The circle graph shows the number of EP-DEGs enriched in Kyoto Encyclopedia of Genes and Genomes (KEGG). The hsa number is the ID value of the KEGG entry. Each column in the inner circle corresponds to an entry, and the higher the column height, the smaller the p-value. The outer circle is the molecule contained in the entry, and different heights represent the logFC. BP, biological process; CC, cellular component; MF, molecular function.
Table II.
Functional enrichment analysis of extracellular protein-differentially expressed genes.
| Ontology | ID | Description | GeneRatio | BgRatio | p-value* | p.adjust | q-value |
|---|---|---|---|---|---|---|---|
| BP | GO:0030198 | ECM organization | 23/55 | 368/18670 | 3.17e-25 | 3.38e-22 | 2.53e-22 |
| BP | GO:0043062 | Extracellular structure organization | 23/55 | 422/18670 | 7.41e-24 | 3.94e-21 | 2.95e-21 |
| BP | GO:0038065 | Collagen-activated signalling pathway | 3/55 | 13/18670 | 6.78e-06 | 4.81e-04 | 3.59e-04 |
| BP | GO:0007568 | Ageing | 7/55 | 321/18670 | 4.16e-05 | 0.001 | 0.001 |
| BP | GO:0098868 | Bone growth | 2/55 | 47/18670 | 0.008 | 0.062 | 0.046 |
| CC | GO:0062023 | Collagen-containing ECM | 32/56 | 406/19717 | 8.89e-40 | 9.43e-38 | 6.46e-38 |
| CC | GO:0005788 | Endoplasmic reticulum lumen | 21/56 | 309/19717 | 5.20e-24 | 2.75e-22 | 1.89e-22 |
| CC | GO:0005604 | Basement membrane | 13/56 | 95/19717 | 5.13e-19 | 1.81e-17 | 1.24e-17 |
| CC | GO:0043202 | Lysosomal lumen | 3/56 | 95/19717 | 0.002 | 0.014 | 0.010 |
| CC | GO:0005796 | Golgi lumen | 3/56 | 102/19717 | 0.003 | 0.016 | 0.011 |
| MF | GO:0005201 | ECM structural constituent | 24/54 | 163/17697 | 2.66e-35 | 3.40e-33 | 2.54e-33 |
| MF | GO:0005518 | Collagen binding | 11/54 | 67/17697 | 8.15e-17 | 5.22e-15 | 3.90e-15 |
| MF | GO:0048407 | Platelet-derived growth factor binding | 5/54 | 11/17697 | 9.97e-11 | 3.19e-09 | 2.39e-09 |
| MF | GO:0019838 | Growth factor binding | 8/54 | 137/17697 | 8.10e-09 | 2.07e-07 | 1.55e-07 |
| MF | GO:0048306 | Calcium-dependent protein binding | 2/54 | 61/17697 | 0.015 | 0.073 | 0.055 |
-
*
Independent-samples t-test.
-
BP, biological process; CC, cellular component; ECM, extracellular matrix; GO, Gene Ontology; MF, molecular function.
In addition, KEGG pathway enrichment analysis on the EP-DEGs showed that the PI3K/Akt signalling pathway, focal adhesion, ECM-receptor interaction, protein digestion and absorption, the AGE/RAGE signalling pathway, the relaxin signalling pathway, platelet activation, and regulation of the actin cytoskeleton were enriched (Supplementary Figure ab and Table III). The circle graph shows the genes enriched in each pathway (Figure 4d).
Table III.
Functional enrichment analysis of extracellular protein-differentially expressed genes.
| Ontology | ID | Description | GeneRatio | BgRatio | p-value* | p.adjust | q-value |
|---|---|---|---|---|---|---|---|
| KEGG | hsa04512 | ECM-receptor interaction | 10/24 | 88/8076 | 2.42e-14 | 1.04e-12 | 7.38e-13 |
| KEGG | hsa04510 | Focal adhesion | 11/24 | 201/8076 | 3.25e-12 | 6.99e-11 | 4.96e-11 |
| KEGG | hsa04974 | Protein digestion and absorption | 9/24 | 103/8076 | 6.99e-12 | 1.00e-10 | 7.11e-11 |
| KEGG | hsa04151 | PI3K-Akt signalling pathway | 11/24 | 354/8076 | 1.47e-09 | 1.58e-08 | 1.12e-08 |
| KEGG | hsa04933 | AGE-RAGE signalling pathway in diabetic complications | 6/24 | 100/8076 | 3.48e-07 | 2.14e-06 | 1.52e-06 |
| KEGG | hsa04926 | Relaxin signalling pathway | 5/24 | 129/8076 | 3.20e-05 | 1.72e-04 | 1.22e-04 |
| KEGG | hsa04611 | Platelet activation | 4/24 | 124/8076 | 4.43e-04 | 0.002 | 0.001 |
| KEGG | hsa04810 | Regulation of actin cytoskeleton | 3/24 | 218/8076 | 0.026 | 0.079 | 0.056 |
-
*
Independent-samples t-test.
-
ECM, extracellular matrix; KEGG, Kyoto Encyclopedia of Genes and Genomes.
Protein-protein interaction network construction and identification of hub EP-DEGs
Using the STRING database, we constructed a PPI network of 56 EP-DEGs and visualized them to investigate the interactions between the corresponding EP-DEG proteins (Supplementary Figure b). The PPI network consists of 56 nodes and 225 edges. Two significant modules (score ≥ 4.5) were extracted from the PPI network. Module 1 contained 19 upregulated gene nodes and 128 edges. Module 2 contained 12 upregulated gene nodes and 34 edges. The darker the colour and the wider the edge, the more substantial the evidence of PPIs (Figures 5a and 5b). We used the cytoHubba in Cytoscape to screen the top ten hub EP-DEGs (Figure 6a, Table IV). The hub EP-DEG expression levels of COL1A1, COL1A2, COL3A1, COL4A2, COL4A1, COL6A3, LUM, FN1, and COL15A1 were markedly upregulated in the late IDD group (Figure 6b). In addition, we performed Pearson’s correlation analysis to explore correlations between the expression levels of the ten hub genes. There were significant positive correlations between all the hub genes except for LUM, COL4A1, and COL4A2, which were not statistically significant (Figure 6c).
Fig. 5
The module analysis of the protein-protein interaction (PPI) network. a) Module 1 contained 19 gene nodes and 128 edges, MCODE score = 6.18. b) Module 2 contained 12 upregulated gene nodes and 34 edges, MCODE score = 14.22.
Fig. 6
The top ten hub extracellular protein-differentially expressed genes (EP-DEGs). a) The ten nodes are shown, ranging from red (high degree value) to yellow (low degree value). b) Expression level analysis of the ten hub genes in both groups. c) Bitmap of the Pearson’s correlation analysis between the ten hub genes. Red denotes a positive correlation, green denotes a negative correlation. IDD, intervertebral disc degeneration. All p-values were calculated using independent-samples t-test.
Table IV.
The information of top ten hub extracellular protein-differentially expressed genes.
| Gene name | Gene description | Expression level | Rank |
|---|---|---|---|
| COL1A1 | Collagen type I α 1 chain | upregulated | 1 |
| COL1A2 | Collagen type I α 2 chain | upregulated | 2 |
| COL3A1 | Collagen type III α 1 chain | upregulated | 3 |
| COL4A2 | Collagen type IV α 2 chain | upregulated | 4 |
| COL4A1 | Collagen type IV α 1 chain | upregulated | 5 |
| COL6A3 | Collagen type VI α 3 chain | upregulated | 6 |
| COL2A1 | Collagen type II α 1 chain | upregulated | 7 |
| LUM | Lumican | upregulated | 8 |
| FN1 | Fibronectin 1 | upregulated | 9 |
| COL15A1 | Collagen type XV α 1 chain | upregulated | 10 |
Clinical nucleus pulposus surgical sample collection and quantitative real-time polymerase chain reaction analysis
We performed qRT-PCR analysis on the top ten hub genes to verify further the expression changes in the early and late IDD groups. First, we collected clinical surgical nucleus pulposus tissue specimens due to the limitations of relying solely on the Pfirrmann grading system to assess the degree of IDD under MRI.20,35 We screened six samples of nucleus pulposus tissue from the early and late IDD groups by MRI combined with a histological examination (H&E staining, Safranin O–Fast Green staining, and Toluidine blue staining). Compared with the early IDD group, the late IDD group showed a substantial loss of nucleus pulposus ECM, resulting in decreased water content and low signal on the MRI T2-weighted images. H&E staining showed that the collagen fibres were neatly arranged in the early IDD group, the nucleus pulposus ECM was homogeneously and deeply stained, and a few small vacuolar cells appeared. However, in the late IDD group the collagen fibres were missing, the arrangement was highly disordered, and many nucleus pulposus cells were present, alongside aggregates, which appeared as large vacuolated cells. Safranin O–Fast Green staining was used to identify the presence of proteoglycan matrix and Toluidine blue staining to detect the accumulation of glycosaminoglycans.36 Safranin O–Fast Green and Toluidine blue staining indicated that ECM proteins (proteoglycan matrix and glycosaminoglycans) were primarily lost in the late IDD group. These results guaranteed the accuracy of our sample collection. Next, we performed qRT-PCR analysis to detect expression changes of the hub genes in both groups of samples. Expression of COL1A1, COL1A2, COL3A1, COL4A2, COL4A1, COL6A3, LUM, FN1, and COL15A1 was significantly increased in the late IDD group, while COL2A1 expression was significantly decreased. These results further demonstrate the accuracy of our analysis (Figure 7).
Fig. 7
Clinical surgical nucleus pulposus tissue sample collection and quantitative real-time polymerase chain reaction (qRT-PCR) analysis. a) Sagittal MRI T2-weighted images of the spine in intervertebral disc degeneration (IDD) patients (scale bar = 25 μm). b) Haematoxylin and eosin (H&E) staining of nucleus pulposus tissue of both groups of patients. Magnification: 40× (scale bar = 25 μm). c) Safranin O–Fast Green staining of nucleus pulposus tissue of both groups of patients. Magnification: 40× (scale bar = 25 μm). d) Toluidine blue staining of nucleus pulposus tissue of both groups of patients. Magnification: 40× (scale bar = 25 μm). e) The qRT-PCR analysis of top ten hub genes. Data were presented as the mean and standard deviation (SD). Independent-samples t-test was used for the statistical analysis of the two groups (**p < 0.01, ***p < 0.001).
TF–gene interactions
NetworkAnalyst was used to investigate the interaction between hub EP-DEGs and TFs. TF–gene interactions of hub EP-DEGs were identified, and interactions with hub EP-DEGs are visualized in Figure 8. The network contains a total of 53 TF–gene interactions with 53 nodes and 82 edges. LUM is regulated by 16 TF–genes, FN1 is regulated by 14 TF–genes, and COL15A1 is regulated by 11 TF–genes. All of these TF–genes regulate more than one hub EP-DEG of the network, which indicates high interaction of the TF–genes with hub EP-DEGs (Figure 8).
Fig. 8
Transcription factor (TF)–gene interaction network with hub extracellular protein-differentially expressed genes (EP-DEGs).
TF–miRNA coregulatory network
The TF–miRNA coregulatory network was used to analyze the interactions of TF genes, miRNAs, and hub EP-DEGs, which may explain the changes in expression of the hub EP-DEGs. NetworkAnalyst was used to generate a TF–miRNA coregulatory interaction network. This interaction can be the upstream factor for regulating the expression of the hub EP-DEGs. The TF–miRNA coregulatory network comprises 218 nodes and 334 edges. A total of 140 miRNAs and 53 TF–genes have interacted with the hub EP-DEGs (Supplementary Figure c).
Identification of candidate drugs
The Enrichr platform was used to identify potential drugs targeting ten hub EP-DEGs based on transcriptome signatures in the DSigDB. According to the p-values, the top ten compounds were selected as candidate drugs. The results showed that electrocorundum CTD 00005364, phenytoin CTD 00006527, and sorafenib CTD 00004146 are the three drug molecules that most EP-DEGs interact with. As these signature drugs were detected for the hub EP-DEGs, they represent common drugs for IDD (Supplementary Figure da). Supplementary Figure db shows the ten candidate drugs for hub EP-DEGs obtained from the DSigDB.
Discussion
IDD is a chronic progressive degenerative spinal disease, the incidence of which positively correlates with age.36,37 Chronic progression of IDD is difficult to reverse with currently available treatments. Although the mechanisms underlying IDD development are known, they remain poorly understood.36,38 Genomics and transcriptomics studies have explored potential therapeutic targets for many diseases, including IDD.39 This study used the GSE23130 dataset from the GEO database to identify the DEGs in two patient groups, followed by the HPA and Uniport databases to screen out EP-DEGs. Next, we performed functional enrichment analysis of differential genes, PPI network construction, and identification of hub EP-DEGs. Finally, we used qRT-PCR to validate further the hub EP-DEGs obtained.
We performed a comprehensive analysis of the GEO, HPA, and Uniport databases and screened 56 EP-DEGs, including 54 upregulated and two downregulated genes. GO analysis revealed that EP-DEGs were markedly enriched in extracellular structure organization, collagen-containing ECM, ageing, endoplasmic reticulum lumen, growth factor binding, and platelet-derived growth factor binding. Under normal conditions, intervertebral disc ECM anabolism and catabolism are in dynamic equilibrium, but disruption of this balance often leads to IDD.38 During IDD, proteoglycan and type II collagen content markedly decreases, while type I collagen content markedly increases, resulting in decreased water absorption, mechanical microenvironment dysfunction, and decreased resistance to ECM loading.12 Increasing evidence suggests that multiple proteases (matrix metalloproteinases (MMPs), a disintegrin and metalloproteinase with thrombospondin motifs (ADAMTSs), cathepsin proteases, etc.) are greatly elevated in degenerated IVD tissues, and play a key role in the hydrolysis of ECM proteins.2,12,38 Ageing results in the accumulation of senescent cells in the IVD, and is a major risk factor for IDD progression.40 The number of senescent NP cells increases substantially during IDD, and senescent cells secrete metabolic factors such as proinflammatory cytokines, matrix-degrading proteases, growth factors, and chemokines, leading to changes in the ECM.41-43 IVD tissue has robust ECM protein synthesis and secretion capabilities and a complex microenvironment that undergoes changes in mechanical load, high osmotic pressure, low pH value, and oxidative stress. Therefore, the endoplasmic reticulum in IVD cells is highly susceptible to external stimulation.44 Increased expression of endoplasmic reticulum stress-related proteins and apoptosis markers, such as GRP78, CHOP, and caspase 12, has been reported in NP tissues from IDD patients.45,46
The enriched KEGG pathways included the PI3K/Akt signalling pathway, focal adhesion, ECM-receptor interaction, AGE/RAGE signalling pathway, and regulation of the actin cytoskeleton. The PI3K/Akt signalling pathway is involved in many biological processes such as proliferation, senescence, and apoptosis of IVD cells.47 Activation of the PI3K/Akt pathway delays IDD progression by increasing the expression of SOX9 to promote the expression of ECM-related proteins in NP cells.48-50 Wang et al51 showed that cyclic mechanical stretching ameliorated the degeneration of nucleus pulposus cells via the PI3K/AKT signalling pathway. Oxidative stress and inflammatory responses mediated by the AGE/RAGE signalling pathway are important in IDD progression.52 Illien-Jünger et al53 suggested that activation of the AGE/RAGE pathway induces hypertrophy and osteogenic differentiation of IVD cells. More importantly, the accumulation of AGEs caused changes in the oxidative microenvironment of nucleus pulposus tissues, which in turn led to IDD. Therefore, targeting the AGE/RAGE pathway may be a promising therapeutic target to delay the progression of IDD. Related studies have shown that the cytoskeleton acts as a mechanotransductor between the ECM and IVD cells to sense changes in the local microenvironment of the IVD.54,55 The cytoskeleton in IVD includes F-actin, tubulin microtubules, and vimentin filaments, which are essential for cell motility, division, protein trafficking, and the secretion of cytokines.56 Tessier et al57 showed that tonicity-responsive enhancer-binding protein (TonEBP) deficiency leads to tissue fibrosis associated with changes in the actin cytoskeleton and adhesion molecules. Ke et al55 showed that inhibiting the RhoA/ROCK1 pathway attenuates compressive stress-induced senescence in human nucleus pulposus cells by regulating actomyosin cytoskeleton remodelling.
We performed a PPI network analysis to identify hub EP-DEGs in IDD. We identified ten EP-DEGs, namely COL1A1, COL1A2, COL3A1, COL4A1, COL4A2, COL6A3, COL2A1, LUM, FN1, and COL15A1. We verified the altered expression of EP-DEGs in IDD by qRT-PCR. Healthy nucleus pulposus has a gel-like structure and is mainly composed of type II collagen. Annulus fibrosus is tough, and has fibres composed of collagen types I and II. During IDD, the type I collagen content notably increased in nucleus pulposus tissues, promoting fibrosis in NP tissues.58 COL1A1, COL1A2, COL3A1, and FN1 are widely used as markers of fibrosis.59,60 Since COL1A1 and COL1A2 are involved in regulating type I collagen synthesis, the increased expression of COL1A1 and COL1A2 during IDD increases type I collagen production, resulting in a blurred boundary between annulus fibrosus and nucleus pulposus.61 Using a rabbit IDD model, Lv et al60 showed that the expression level of COL3A1 was greatly increased in degenerated rabbit nucleus pulposus tissue, which was related to the fibrosis of nucleus pulposus tissue. According to Sun et al,62 the expression level of FN1 was markedly increased in degenerated nucleus pulposus cells. LUM is an important ECM glycoprotein in human cartilage tissue and is highly expressed in nucleus pulposus specimens from patients with lumbar disc herniation.63 Li et al64 showed high expression of LUM in IDD, and LUM knockdown effectively attenuated tumour necrosis factor (TNF)-α-induced inflammatory response, cell cycle arrest, and cellular senescence. COL2A1 is a major component of type II collagen and plays an important role in maintaining the gel-like properties of nucleus pulposus tissues. COL2A1 expression is substantially decreased during IDD.65 In this study, the expression level of COL2A1 did not decrease substantially, which may be related to the low sample size. COL4A1 and COL4A2 are the main components of type IV collagen fibres. Type VI collagen is a non-fibrillar collagen, and COL4A1-COL4A2-COL4A3 heterotrimers form bead-like microfibrils that are important for anchoring the basement membrane into the ECM components.66 The COL6A3 gene encodes the α3 chain of type VI collagen fibres and is also an important part of the basement membrane. Dysregulation of the COL6A3 gene can lead to muscular dystrophy.67 Unfortunately, studies of COL4A1, COL4A2, and COL6A3 in IDD are lacking to date.
TF–gene interaction network was obtained with the hub EP-DEGs. TF–genes work as regulators according to genetic expressions, which may result in IDD. Among the regulators, FOXC1, STAT3, FOXL1, and STAT1 have a notable interaction. Many studies have shown that STAT3 and SATA1 upregulate and mediate different biological processes in IDD.68-70 Regulatory biomolecules may act as key factors in IDD. The activities of miRNAs and TF–genes are analyzed for the regulation of hub EP-DEGs that are visualized in the TF–miRNA coregulatory network. Li et al71 showed that ATF3 positively regulates the apoptosis, reactive oxygen species (ROS) production, inflammatory response, and ECM degradation of nucleus pulposus cells induced by tert-butyl hydrogen peroxide. MiR-200c inhibits X-linked inhibitor of apoptosis protein (XiAP)-induced nucleus pulposus cell apoptosis and extrinsic interstitial anabolic and catabolic imbalance.72 Interleukin-1β and TNF-α exposure greatly reduced miR-578 levels in nucleus pulposus, where ectopic miR-578 expression inhibited cell growth, proinflammatory cytokine expression, and ECM degradation.73 Moreover, the DSigDB database was used to predict several drugs, which are related to the hub EP-DEGs. Among all the candidate drugs, this study shows the top ten most notable drugs, which includes electrocorundum, phenytoin, cytarabine, testosterone enanthate, paclitaxel alpha-Neu5Ac, acetaldehyde, and sorafenib. Studies have shown that high concentrations of phenytoin notably attenuate lipopolysaccharide-induced TNF-α production, showing strong anti-ECM catabolic effects.74
Our study has some limitations. First, the bioinformatics analysis results are based on pre-set criteria, and changing the cut-off criteria will ultimately affect the results. Second, due to a lack of IDD datasets of sufficient sample size, it is impossible to validate our analysis results on different datasets. Moreover, although we assessed changes in the hub EP-DEGs at the gene level, alterations in the protein expression levels of these genes were not assessed. We will investigate the role of these hub EP-DEGs in the occurrence and development of IDD in future studies.
In conclusion, we identified COL1A1, COL1A2, COL3A1, COL4A2, COL4A1, COL6A3, LUM, FN1, and COL15A1 as IDD-related hub EP-DEGs through bioinformatics analysis. Our findings provide new insights into IDD pathogenesis and treatment opinions. However, further studies are needed to validate our results.
References
1. Li G , Ma L , He S , et al. WTAP-mediated m6A modification of lncRNA NORAD promotes intervertebral disc degeneration . Nat Commun . 2022 ; 13 ( 1 ): 1469 . Crossref PubMed Google Scholar
2. Zhang G-Z , Chen H-W , Deng Y-J , et al. BRD4 Inhibition Suppresses Senescence and Apoptosis of Nucleus Pulposus Cells by Inducing Autophagy during Intervertebral Disc Degeneration: An In Vitro and In Vivo Study . Oxid Med Cell Longev . 2022 ; 2022 : 9181412 . Crossref PubMed Google Scholar
3. Hu B , Wang P , Zhang S , et al. HSP70 attenuates compression-induced apoptosis of nucleus pulposus cells by suppressing mitochondrial fission via upregulating the expression of SIRT3 . Exp Mol Med . 2022 ; 54 ( 3 ): 309 – 323 . Crossref PubMed Google Scholar
4. Zhang Z , Wang Q , Li Y , Li B , Zheng L , He C . Hub Genes and Key Pathways of Intervertebral Disc Degeneration: Bioinformatics Analysis and Validation . Biomed Res Int . 2021 ; 2021 : 5340449 . Crossref PubMed Google Scholar
5. Ohnishi T , Iwasaki N , Sudo H . Causes of and Molecular Targets for the Treatment of Intervertebral Disc Degeneration: A Review . Cells . 2022 ; 11 ( 3 ): 394 . Crossref PubMed Google Scholar
6. Kuyl E-V , Shu F , Sosa BR , et al. Inhibition of PAD4 mediated neutrophil extracellular traps prevents fibrotic osseointegration failure in a tibial implant murine model: an animal study . Bone Joint J . 2021 ; 103-B ( 7 Supple B ): 135 – 144 . Crossref PubMed Google Scholar
7. Chen H-W , Zhang G-Z , Liu M-Q , et al. Natural Products of Pharmacology and Mechanisms in Nucleus Pulposus Cells and Intervertebral Disc Degeneration . Evid Based Complement Alternat Med . 2021 ; 2021 : 9963677 . Crossref PubMed Google Scholar
8. Wu Z-L , Xie Q-Q , Liu T-C , Yang X , Zhang G-Z , Zhang H-H . Role of the Wnt pathway in the formation, development, and degeneration of intervertebral discs . Pathol Res Pract . 2021 ; 220 : 153366 . Crossref PubMed Google Scholar
9. Roh EJ , Darai A , Kyung JW , et al. Genetic Therapy for Intervertebral Disc Degeneration . Int J Mol Sci . 2021 ; 22 ( 4 ): 1579 . Crossref PubMed Google Scholar
10. Yang X , Arts MP , Bartels R , Vleggeert-Lankamp CLA . The type of cervical disc herniation on MRI does not correlate to clinical outcomes . Bone Joint J . 2022 ; 104-B ( 11 ): 1242 – 1248 . Crossref PubMed Google Scholar
11. Lang G , Liu Y , Geries J , et al. An intervertebral disc whole organ culture system to investigate proinflammatory and degenerative disc disease condition . J Tissue Eng Regen Med . 2018 ; 12 ( 4 ): e2051 – e2061 . Crossref PubMed Google Scholar
12. Liang H , Luo R , Li G , Zhang W , Song Y , Yang C . The Proteolysis of ECM in Intervertebral Disc Degeneration . Int J Mol Sci . 2022 ; 23 ( 3 ): 1715 . Crossref PubMed Google Scholar
13. Liu Z-M , Lu C-C , Shen P-C , et al. Suramin attenuates intervertebral disc degeneration by inhibiting NF-κB signalling pathway . Bone Joint Res . 2021 ; 10 ( 8 ): 498 – 513 . Crossref PubMed Google Scholar
14. Francisco V , Pino J , González-Gay MÁ , et al. A new immunometabolic perspective of intervertebral disc degeneration . Nat Rev Rheumatol . 2022 ; 18 ( 1 ): 47 – 60 . Crossref PubMed Google Scholar
15. Lai MKL , Cheung PWH , Samartzis D , Karppinen J , Cheung KMC , Cheung JPY . The profile of the spinal column in subjects with lumbar developmental spinal stenosis . Bone Joint J . 2021 ; 103-B ( 4 ): 725 – 733 . Crossref PubMed Google Scholar
16. Liu Z-M , Shen P-C , Lu C-C , Chou S-H , Tien Y-C . Suramin enhances chondrogenic properties by regulating the p67phox/PI3K/AKT/SOX9 signalling pathway . Bone Joint Res . 2022 ; 11 ( 10 ): 723 – 738 . Crossref PubMed Google Scholar
17. Edgar R , Domrachev M , Lash AE . Gene Expression Omnibus: NCBI gene expression and hybridization array data repository . Nucleic Acids Res . 2002 ; 30 ( 1 ): 207 – 210 . Crossref PubMed Google Scholar
18. No authors listed . Gene Expression Omnibus (GEO). National Center for Biotechnology Information (NCBI) . 2023 . https://www.ncbi.nlm.nih.gov/geo/ ( date last accessed 11 July 2023 ). Google Scholar
19. Davis S , Meltzer PS . GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor . Bioinformatics . 2007 ; 23 ( 14 ): 1846 – 1847 . Crossref PubMed Google Scholar
20. Pfirrmann CW , Metzdorf A , Zanetti M , Hodler J , Boos N . Magnetic resonance classification of lumbar intervertebral disc degeneration . Spine (Phila Pa 1976) . 2001 ; 26 ( 17 ): 1873 – 1878 . Crossref PubMed Google Scholar
21. No authors listed . Bioconductor . R programming language . 2003 . https://www.bioconductor.org/ ( date last accessed 11 July 2023 ). Google Scholar
22. Smyth GK . Limma: linear models for microarray data . In : Gentleman R , Carey V , Dudoit S , Irizarry R , Huber W . Bioinformatics and Computational Biology Solutions Using R and Bioconductor . New York, New York : Springer , 2005 : 397 – 420 . Crossref Google Scholar
23. No authors listed . The Human Protein Atlas . 2023 . https://www.proteinatlas.org/ ( date last accessed 11 July 2023 ). Crossref PubMed Google Scholar
24. No authors listed . Uniprot . UniProt consortium . 2002 . https://www.uniprot.org/ ( date last accessed 11 July 2023 ). Crossref PubMed Google Scholar
25. Uhlén M , Fagerberg L , Hallström BM , et al. Proteomics. Tissue-based map of the human proteome . Science . 2015 ; 347 ( 6220 ): 1260419 . Crossref PubMed Google Scholar
26. UniProt Consortium . UniProt: the universal protein knowledgebase in 2021 . Nucleic Acids Res . 2021 ; 49 ( D1 ): D480 – D489 . Crossref PubMed Google Scholar
27. Yu G , Wang L-G , Han Y , He Q-Y . clusterProfiler: an R package for comparing biological themes among gene clusters . OMICS . 2012 ; 16 ( 5 ): 284 – 287 . Crossref PubMed Google Scholar
28. No authors listed . STRING: Search Tool for the Retrieval of Interacting Genes/Proteins . STRING Consortium . 2023 . https://cn.string-db.org/ ( date last accessed 11 July 2023 ). Google Scholar
29. Szklarczyk D , Gable AL , Lyon D , et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets . Nucleic Acids Res . 2019 ; 47 ( D1 ): D607 – D613 . Crossref PubMed Google Scholar
30. No authors listed . Cytoscape . Cytoscape Consortium . 2001 . https://cytoscape.org/ ( date last accessed 11 July 2023 ). Crossref PubMed Google Scholar
31. Chin C-H , Chen S-H , Wu H-H , Ho C-W , Ko M-T , Lin C-Y . cytoHubba: identifying hub objects and sub-networks from complex interactome . BMC Syst Biol . 2014 ; 8 Suppl 4 ( Suppl 4 ): S11 . Crossref PubMed Google Scholar
32. No authors listed . NetworkAnalyst . 2023 . https://www.networkanalyst.ca/ ( date last accessed 11 July 2023 ). Crossref PubMed Google Scholar
33. No authors listed . JASPAR - A database of transcription factor binding profiles . 2023 . https://jaspar.genereg.net/ ( date last accessed 11 July 2023 ). Google Scholar
34. No authors listed . Enrichr . 2023 . https://maayanlab.cloud/Enrichr/ ( date last accessed 19 July 2023 ). Crossref PubMed Google Scholar
35. Griffith JF , Wang Y-X , Antonio GE , et al. Modified Pfirrmann grading system for lumbar intervertebral disc degeneration . Spine . 2007 ; 32 ( 24 ): E708 – 12 . Crossref PubMed Google Scholar
36. Zhang G-Z , Deng Y-J , Xie Q-Q , et al. Sirtuins and intervertebral disc degeneration: Roles in inflammation, oxidative stress, and mitochondrial function . Clin Chim Acta . 2020 ; 508 : 33 – 42 . Crossref PubMed Google Scholar
37. Zhu D , Zhou W , Wang Z , et al. Periostin: An Emerging Molecule With a Potential Role in Spinal Degenerative Diseases . Front Med (Lausanne) . 2021 ; 8 : 694800 . Crossref PubMed Google Scholar
38. Zhang G-Z , Liu M-Q , Chen H-W , et al. NF-κB signalling pathways in nucleus pulposus cell function and intervertebral disc degeneration . Cell Prolif . 2021 ; 54 ( 7 ): e13057 . Crossref PubMed Google Scholar
39. Wang X , Tan J , Sun J , et al. Transcriptomics Study to Determine the Molecular Mechanism by which sIL-13Rα2-Fc Inhibits Caudal Intervertebral Disc Degeneration in Rats . Biomed Res Int . 2020 ; 2020 : 7645989 . Crossref PubMed Google Scholar
40. Wang F , Cai F , Shi R , Wang X-H , Wu X-T . Aging and age related stresses: a senescence mechanism of intervertebral disc degeneration . Osteoarthr Cartil . 2016 ; 24 ( 3 ): 398 – 408 . Crossref PubMed Google Scholar
41. Dudek M , Yang N , Ruckshanthi JP , et al. The intervertebral disc contains intrinsic circadian clocks that are regulated by age and cytokines and linked to degeneration . Ann Rheum Dis . 2017 ; 76 ( 3 ): 576 – 584 . Crossref PubMed Google Scholar
42. Sakai D , Nakamura Y , Nakai T , et al. Exhaustion of nucleus pulposus progenitor cells with ageing and degeneration of the intervertebral disc . Nat Commun . 2012 ; 3 : 1264 . Crossref PubMed Google Scholar
43. Wu Y , Shen S , Shi Y , Tian N , Zhou Y , Zhang X . Senolytics: Eliminating Senescent Cells and Alleviating Intervertebral Disc Degeneration . Front Bioeng Biotechnol . 2022 ; 10 : 823945 . Crossref PubMed Google Scholar
44. Wang D , He X , Zheng C , et al. Endoplasmic Reticulum Stress: An Emerging Therapeutic Target for Intervertebral Disc Degeneration . Front Cell Dev Biol . 2021 ; 9 : 819139 . Crossref PubMed Google Scholar
45. Liao Z , Luo R , Li G , et al. Exosomes from mesenchymal stem cells modulate endoplasmic reticulum stress to protect against nucleus pulposus cell death and ameliorate intervertebral disc degeneration in vivo . Theranostics . 2019 ; 9 ( 14 ): 4084 – 4100 . Crossref PubMed Google Scholar
46. Zhao C-Q , Zhang Y-H , Jiang S-D , Jiang L-S , Dai L-Y . Both endoplasmic reticulum and mitochondria are involved in disc cell apoptosis and intervertebral disc degeneration in rats . Age (Dordr) . 2010 ; 32 ( 2 ): 161 – 177 . Crossref PubMed Google Scholar
47. Li J , Wang C , Xue L , Zhang F , Liu J . Melatonin Suppresses Apoptosis of Nucleus Pulposus Cells through Inhibiting Autophagy via the PI3K/Akt Pathway in a High-Glucose Culture . Biomed Res Int . 2021 ; 2021 : 4604258 . Crossref PubMed Google Scholar
48. Tan Y , Yao X , Dai Z , Wang Y , Lv G . Bone morphogenetic protein 2 alleviated intervertebral disc degeneration through mediating the degradation of ECM and apoptosis of nucleus pulposus cells via the PI3K/Akt pathway . Int J Mol Med . 2019 ; 43 ( 1 ): 583 – 592 . Crossref PubMed Google Scholar
49. Cheng C-C , Uchiyama Y , Hiyama A , Gajghate S , Shapiro IM , Risbud MV . PI3K/AKT regulates aggrecan gene expression by modulating Sox9 expression and activity in nucleus pulposus cells of the intervertebral disc . J Cell Physiol . 2009 ; 221 ( 3 ): 668 – 676 . Crossref PubMed Google Scholar
50. Zhang X , Hu Z , Hao J , Shen J . Low Intensity Pulsed Ultrasound Promotes the Extracellular Matrix Synthesis of Degenerative Human Nucleus Pulposus Cells Through FAK/PI3K/Akt Pathway . Spine . 2016 ; 41 ( 5 ): E248 – 54 . Crossref PubMed Google Scholar
51. Wang D , Chen Y , Cao S , et al. Cyclic Mechanical Stretch Ameliorates the Degeneration of Nucleus Pulposus Cells through Promoting the ITGA2/PI3K/AKT Signaling Pathway . Oxid Med Cell Longev . 2021 ; 2021 : 6699326 . Crossref PubMed Google Scholar
52. Lin J , Gu J , Fan D , Li W . Herbal Formula Modified Bu-Shen-Huo-Xue Decoction Attenuates Intervertebral Disc Degeneration via Regulating Inflammation and Oxidative Stress . Evid Based Complement Alternat Med . 2022 ; 2022 : 4284893 . Crossref PubMed Google Scholar
53. Illien-Jünger S , Torre OM , Kindschuh WF , Chen X , Laudier DM , Iatridis JC . AGEs induce ectopic endochondral ossification in intervertebral discs . Eur Cell Mater . 2016 ; 32 : 257 – 270 . Crossref PubMed Google Scholar
54. Tsingas M , Ottone OK , Haseeb A , et al. Sox9 deletion causes severe intervertebral disc degeneration characterized by apoptosis, matrix remodeling, and compartment-specific transcriptomic changes . Matrix Biol . 2020 ; 94 : 110 – 133 . Crossref PubMed Google Scholar
55. Ke W , Wang B , Hua W , et al. The distinct roles of myosin IIA and IIB under compression stress in nucleus pulposus cells . Cell Prolif . 2021 ; 54 ( 2 ): e12987 . Crossref PubMed Google Scholar
56. Zhang C , Wang F , Xie Z , Chen L , Wu X . The hippo pathway orchestrates cytoskeletal organisation during intervertebral disc degeneration . Acta Histochem . 2021 ; 123 ( 6 ): 151770 . Crossref PubMed Google Scholar
57. Tessier S , Tran VA , Ottone OK , et al. TonEBP-deficiency accelerates intervertebral disc degeneration underscored by matrix remodeling, cytoskeletal rearrangements, and changes in proinflammatory gene expression . Matrix Biol . 2020 ; 87 : 94 – 111 . Crossref PubMed Google Scholar
58. Gao X , Zhu Q , Gu W . Prediction of glycosaminoglycan synthesis in intervertebral disc under mechanical loading . J Biomech . 2016 ; 49 ( 13 ): 2655 – 2661 . Crossref PubMed Google Scholar
59. Yang X , Chen Z , Chen C , et al. Bleomycin induces fibrotic transformation of bone marrow stromal cells to treat height loss of intervertebral disc through the TGFβR1/Smad2/3 pathway . Stem Cell Res Ther . 2021 ; 12 ( 1 ): 34 . Crossref PubMed Google Scholar
60. Lv F-J , Peng Y , Lim FL , et al. Matrix metalloproteinase 12 is an indicator of intervertebral disc degeneration co-expressed with fibrotic markers . Osteoarthr Cartil . 2016 ; 24 ( 10 ): 1826 – 1836 . Crossref PubMed Google Scholar
61. Shen L , Xiao Y , Wu Q , Liu L , Zhang C , Pan X . TLR4/NF-κB axis signaling pathway-dependent up-regulation of miR-625-5p contributes to human intervertebral disc degeneration by targeting COL1A1 . Am J Transl Res . 2019 ; 11 ( 3 ): 1374 – 1388 . PubMed Google Scholar
62. Sun Y , Lv M , Zhou L , et al. Enrichment of committed human nucleus pulposus cells expressing chondroitin sulfate proteoglycans under alginate encapsulation . Osteoarthr Cartil . 2015 ; 23 ( 7 ): 1194 – 1203 . Crossref PubMed Google Scholar
63. Hu S , Fu Y , Yan B , Shen Z , Lan T . Analysis of key genes and pathways associated with the pathogenesis of intervertebral disc degeneration . J Orthop Surg Res . 2020 ; 15 ( 1 ): 371 . Crossref PubMed Google Scholar
64. Li Z , Sun C , Chen M , Wang B . Lumican silencing alleviates tumor necrosis factor-α-induced nucleus pulposus cell inflammation and senescence by inhibiting apoptosis signal regulating kinase 1/p38 signaling pathway via inactivating Fas ligand expression . Bioengineered . 2021 ; 12 ( 1 ): 6891 – 6901 . Crossref PubMed Google Scholar
65. Ni L , Zheng Y , Gong T , et al. Proinflammatory macrophages promote degenerative phenotypes in rat nucleus pulpous cells partly through ERK and JNK signaling . J Cell Physiol . 2019 ; 234 ( 5 ): 5362 – 5371 . Crossref PubMed Google Scholar
66. Steffensen LB , Rasmussen LM . A role for collagen type IV in cardiovascular disease? Am J Physiol Heart Circ Physiol . 2018 ; 315 ( 3 ): H610 – H625 . Crossref PubMed Google Scholar
67. Guadagnin E , Mohassel P , Johnson KR , et al. Transcriptome analysis of collagen VI-related muscular dystrophy muscle biopsies . Ann Clin Transl Neurol . 2021 ; 8 ( 11 ): 2184 – 2198 . Crossref PubMed Google Scholar
68. Bai X , Jiang M , Wang J , et al. Cyanidin attenuates the apoptosis of rat nucleus pulposus cells and the degeneration of intervertebral disc via the JAK2/STAT3 signal pathway in vitro and in vivo . Pharm Biol . 2022 ; 60 ( 1 ): 427 – 436 . Crossref PubMed Google Scholar
69. Zhou T , Lin H , Cheng Z , Ji C , Zhang C , Tian J . Mechanism of microRNA-146a-mediated IL-6/STAT3 signaling in lumbar intervertebral disc degeneration . Exp Ther Med . 2017 ; 14 ( 2 ): 1131 – 1135 . Crossref PubMed Google Scholar
70. Hu B , Wang J , Wu X , Chen Y , Yuan W , Chen H . Interleukin-17 upregulates vascular endothelial growth factor by activating the JAK/STAT pathway in nucleus pulposus cells . Joint Bone Spine . 2017 ; 84 ( 3 ): 327 – 334 . Crossref PubMed Google Scholar
71. Li Y , Pan D , Wang X , et al. Silencing ATF3 Might Delay TBHP-Induced Intervertebral Disc Degeneration by Repressing NPC Ferroptosis, Apoptosis, and ECM Degradation . Oxid Med Cell Longev . 2022 ; 2022 : 4235126 . Crossref PubMed Google Scholar
72. Cheng X , Zhang L , Zhang K , et al. Circular RNA VMA21 protects against intervertebral disc degeneration through targeting miR-200c and X linked inhibitor-of-apoptosis protein . Ann Rheum Dis . 2018 ; 77 ( 5 ): 770 – 779 . Crossref PubMed Google Scholar
73. Yan P , Sun C , Luan L , et al. Hsa_circ_0134111 promotes intervertebral disc degeneration via sponging miR-578 . Cell Death Discov . 2022 ; 8 ( 1 ): 55 . Crossref PubMed Google Scholar
74. Serra R , Al-Saidi A-G , Angelov N , Nares S . Suppression of LPS-induced matrix-metalloproteinase responses in macrophages exposed to phenytoin and its metabolite, 5-(p-hydroxyphenyl-), 5-phenylhydantoin . J Inflamm . 2010 ; 7 : 48 . Crossref PubMed Google Scholar
Author contributions
G-Z. Zhang: Investigation, Methodology, Formal analysis, Writing – original draft, Writing – review & editing.
L. Li: Data curation, Formal analysis, Writing – original draft.
Z-B. Luo: Investigation, Formal analysis, Writing – original draft.
C-Y. Zhang: Investigation, Formal analysis, Writing – original draft.
Y-G Wang: Investigation, Formal analysis, Writing – original draft.
X-W. Kang: Supervision, Formal analysis, Writing – review & editing.
Funding statement
The authors disclose receipt of the following financial or material support for the research, authorship, and/or publication of this article: National Natural Sciences Foundation of China (Grant No. 82272536), Natural Science Foundation of Gansu province (Grant No. 21JR7RA406, 22JR5RA949, and 23JRRA1012), the Special Fund Project for Doctoral Training Program of Lanzhou University Second Hospital (Grant No. YJS-BD-09), and Cuiying Scientific and Technological Innovation Program of Lanzhou University Second Hospital (Grant No.CY2020-MS20, and CY2022-QN-A06). Lanzhou University Innovation and Entrepreneurship Cultivation Project (Grant No. cxcy2023012).
ICMJE COI statement
The authors have no conflicts of interest to declare.
Data sharing
The data that support the findings for this study are available to other researchers from the corresponding author upon reasonable request.
Acknowledgements
This study was a re-analysis based on published data from the Gene Expression Omnibus (GEO) database. The authors acknowledge the GEO database for the sharing of data, the Xiantao tool online database (https://www.xiantao.love/) for data visualization, and Textcheck (www.textcheck.com) for the linguistic editing and proofreading services during the preparation of this manuscript.
Ethical review statement
The 12 nucleus pulposus tissue specimens used in this study were collected between June 2019 and December 2020 from patients who underwent spinal surgery in the Orthopedics Department of Lanzhou University Second Hospital, due to intervertebral disc (IVD) herniation or idiopathic scoliosis. This study was performed in compliance with the Declaration of Helsinki and was approved by the Medical Ethics Committee of Lanzhou University Second Hospital (2022A-138).
Open access funding
The authors report that the open access funding for their manuscript was self-funded.
Supplementary material
Figures showing functional enrichment analysis of extracellular protein-differentially expressed genes (EP-DEGs), protein-protein interaction network of EP-DEGs, transcription factor–microRNA coregulatory network with hub EP-DEGs, and candidate drugs for hub EP-DEGs from the Drug Signatures Database. Table showing the details of the 54 upregulated and two downregulated EP-DEGs.
© 2023 Author(s) et al. This is an open-access article distributed under the terms of the Creative Commons Attribution Non-Commercial No Derivatives (CC BY-NC-ND 4.0) licence, which permits the copying and redistribution of the work only, and provided the original author and source are credited. See https://creativecommons.org/licenses/by-nc-nd/4.0/
