An Oxidative Stress-Related Gene Pair (CCNB1/ PKD1), Competitive Endogenous RNAs, and Immune-Infiltration Patterns Potentially Regulate Intervertebral Disc Degeneration Development

Abstract

Oxidative stress (OS) irreversibly affects the pathogenesis of intervertebral disc degeneration (IDD). Certain non-coding RNAs act as competitive endogenous RNAs (ceRNAs) that regulate IDD progression. Analyzing the signatures of oxidative stress-related gene (OSRG) pairs and regulatory ceRNA mechanisms and immune-infiltration patterns associated with IDD may enable researchers to distinguish IDD and reveal the underlying mechanisms. In this study, OSRGs were downloaded and identified using the Gene Expression Omnibus database. Functional-enrichment analysis revealed the involvement of oxidative stress-related pathways and processes, and a ceRNA network was generated. Differentially expressed oxidative stress-related genes (De-OSRGs) were used to construct De-OSRG pairs, which were screened, and candidate De-OSRG pairs were identified. Immune cell-related gene pairs were selected via immune-infiltration analysis. A potential long non-coding RNA-microRNA-mRNA axis was determined, and clinical values were assessed. Eighteen De-OSRGs were identified that were primarily related to intricate signal-transduction pathways, apoptosis-related biological processes, and multiple kinase-related molecular functions. A ceRNA network consisting of 653 long non-coding RNA-microRNA links and 42 mRNA-miRNA links was constructed. Three candidate De-OSRG pairs were screened out from 13 De-OSRG pairs. The abundances of resting memory CD4⁺ T cells, resting dendritic cells, and CD8⁺ T cells differed between the control and IDD groups. CD8⁺ T cell infiltration correlated negatively with cyclin B1 (CCNB1) expression and positively with protein kinase D1 (PKD1) expression. CCNB1-PKD1 was the only pair that was differentially expressed in IDD, was correlated with CD8⁺ T cells, and displayed better predictive accuracy compared to individual genes. The PKD1-miR-20b-5p-AP000797 and CCNB1-miR-212-3p-AC079834 axes may regulate IDD. Our findings indicate that the OSRG pair CCNB1-PKD1, which regulates oxidative stress during IDD development, is a robust signature for identifying IDD. This OSRG pair and increased infiltration of CD8⁺ T cells, which play important roles in IDD, were functionally associated. Thus, the OSRG pair CCNB1-PKD1 is promising target for treating IDD.

Keywords: competitive endogenous RNA; cyclin B1; differential expression; gene pair; immune infiltration; intervertebral disc degeneration; oxidative stress; protein kinase D1.

Figure 1

Flowchart. Triangle represents target‐gene prediction (first triangle, lncRNA target prediction; second triangle, miRNA target prediction).

Figure 2

Volcano plots of DEGs (A–C). (A) Volcano plot of De–lncRNAs. (B) Volcano plot of De-miRNAs. (C) Volcano plot of De-mRNAs. Orange represents up-regulated differentially expressed genes, while purple represents down-regulated differentially expressed genes. The ceRNA network analyses of 18 DEOSRGs, and their function annotation by GO and KEGG (D–G). (D) The ceRNA network analysis graph. GO enrichment analysis of 18 DEOSRGs (E, F). (E) BP, biological process; (F) MF, molecular function. (G) KEGG pathway enrichment analysis of 18 DEOSRGs. De-lncRNAs, differentially expressed long noncoding RNAs; De–miRNAs, differentially expressed microRNAs; De-mRNAs, differentially expressed mRNAs; De–OSRGs, differentially expressed oxidative stress-related genes. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 3

The expression profile for De-OSRGs (A) and De-OSRGPs (B). The difference between OSRGs heat map and OSRGPs heat map (A, B): although the batch correction is used, the gap between different datasets can be clearly shown in panel (A) The screening process of key OSRGPs for IDD. (C) LASSO regression to screen for candidate signatures. (D) SVM to screen for candidate signatures. (E) Venn diagram demonstrating the intersection of candidate signatures screened by LASSO regression and SVM. OSRGs, oxidative stress-related genes; OSRGPs, oxidative stress-related gene pairs; IDD, intervertebral disc degeneration; LASSO, least absolute shrinkage and selection operator; SVM, support vector machine.

Figure 4

Visualization of immune cell infiltration analysis. (A) Violin diagram of the abundance of infiltration by 22 immune cell subsets between IDD and control groups. Blue and orange colors represent IDD and control groups, respectively. (B) Histogram of the distribution of infiltrating immune cells in each individual. Columns represent IDD individuals. (C) Correlation heatmap of 22 immune cell subsets. Blue represents positive correlation, whereas purple represents negative correlation. The color and number in each circle indicate the correlation index between the two immune cell subsets. IDD, intervertebral disc degeneration.

Figure 5

Correlation analysis between infiltrating immune cell subtypes with the candidate De-OSRGs. Spearman correlation analysis between immune cell subtypes and CCNB1 (A) /PKD1 (B) /COX15 (C). Venn diagram (D) showing the key gene pair shared by the immune cell-related genes (ICRG) and candidate De-OSRGs.

Figure 6

|The 69lncRNAs–2miRNAs–2mRNAs (CCNB1/PKD1 pair) ceRNA co-expression network was constructed (A) and was used to screen the key ceRNAs (lncRNAs) by SVM (B, D) and Lasso (C, E). Panels (B, C), venn result for SVM (B) and LASSO (C) is AC079834. LncRNA AC079834 was selected in SVM and LASSO, and was used for building axis (CCNB1–miR-212-3p–AC079834). Panels (D, E), venn result for SVM (D) and LASSO (E) is AP000797. LncRNA AP000797 was selected in SVM and LASSO, and was used for building axis (PKD1–miR-20b-5p–AP000797). Parameters and results in detail are shown in Supplementary Tables S9, S10 .The expression profile for CCNB1–miR-212-3p–AC079834 axis (F) and PKD1–miR-20b-5p–AP000797 axis (G). The low expression of CCNB1 is negatively correlated with increased CD8⁺ T cells infiltration; the high expression of PKD1 is positively correlated with increased CD8⁺ T cells infiltration (H). * represents p < 0.05.

Figure 7

Comparison of ROC analysis between gene pair and separate gene in predicting IDD in the expanded dataset (A–C) and external validation dataset (D–F). ROCs for CCNB1 (A, D), PKD1 (B, E), and CCNB1–PKD1 pairs (C, F).