ABSTRACT Background Low back pain is a significant burden worldwide, and intervertebral disc degeneration (IVDD) is identified as the primary cause. Recent research has emphasized the significant role of circadian rhythms (CRs) and immunity in affecting intervertebral discs (IVD). However, the influence of circadian rhythms and immunity on the mechanism of IVDD remains unclear. This study aimed to identify and validate key rhythm‐related genes in IVDD and analyze their correlation with immune cell infiltration. Methods Two gene expression profiles related to IVDD and rhythm‐related genes were obtained from the Gene Expression Omnibus and GeneCards databases to identify differentially expressed rhythm‐related genes (DERGs). Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA) were conducted to explore the biological functions of these genes. LASSO regression and SVM algorithms were employed to identify hub genes. We subsequently investigated the correlation between hub rhythm‐related genes and immune cell infiltration. Finally, nucleus pulposus‐derived mesenchymal stem cells (NPMSCs) were isolated from normal and degenerative human IVD tissues. Hub rhythm‐related genes expression in NPMSCs was confirmed by real‐time quantitative PCR (RT‐qPCR). Results Six hub genes related to CRs (CCND1, FOXO1, FRMD8, NTRK2, PRRT1, and TFPI) were screened out. Immune infiltration analysis revealed that the IVDD group had significantly more M0 macrophages and significantly fewer follicular helper T cells than those of the control group. Specifically, M0 macrophages were significantly associated with FRMD8, PRRT1, and TFPI. T follicular helper cells were significantly associated with FRDM8, FOXO1, and CCND1. We further confirmed that CCND1, FRMD8, NTRK2, and TFPI were dysrhythmic within NPMSCs from degenerated IVD in vitro. Conclusion Six genes (CCND1, FOXO1, FRMD8, NTRK2, PRRT1 and TFPI) linked to circadian rhythms associated with IVDD progression, together with immunity. The identification of these DEGs may provide new insights for the diagnosis and treatment of IVDD. Keywords: bioinformatics, circadian rhythms, immune cell infiltration, intervertebral disc degeneration, machine learning algorithms, mesenchymal stem cell __________________________________________________________________ Within this manuscript, bioinformatics and experimental methods are used to explore the role of circadian rhythm and immunity in Intervertebral Disc Degeneration (IVDD). CCND1, FOXO1, FRMD8, NTRK2, PRRT1, and TFPI were successfully demonstrated to be hub rhythm‐related genes in the process of IVDD together with immunity. Nucleus pulposus‐derived mesenchymal stem cells (NPMSCs) were used to verify hub genes. graphic file with name JSP2-8-e70066-g008.jpg __________________________________________________________________ Abbreviations AF annulus fibrosus BMAL1 brain and muscle aryl hydrocarbon receptor nuclear translocator‐like 1 CCL2 inflammatory monocyte chemokine ligand CCND1 cyclin D1 CIBERSORT cell‐type identification by estimating relative subset proportions CLOCK circadian locomotor output cycles kaput CRY cryptochrome homologs DEGs differentially expressed genes DERGs differentially expressed rhythm‐related genes FOXO1 forkhead box O1 FRMD8 ferm domain containing 8 GO gene ontology GSEA gene set enrichment analysis HIF‐1α hypoxia inducible factor 1 subunit alpha IL‐6 interleukin‐6 IVDD intervertebral disc degeneration KEGG kyoto encyclopedia of genes and genomes LASSO least absolute shrinkage and selection operator LBP low back pain MHC major histocompatibility complex NP nucleus pulposus NPMSCs nucleus pulposus mesenchymal stem cells NPMSCs nucleus pulposus‐derived mesenchymal stem cells NTRK2 neurotrophic tropomyosin receptor kinase B PBS phosphate‐buffered saline PCA principal component analysis PER period homologs PPAR peroxisome proliferator‐activated receptor PRRT1 proline‐rich transmembrane protein 1 REV‐ERB nuclear receptor subfamily 1 group D ROC receiver operating characteristic ROR retinoid‐related orphan receptor SIRT1 NAD‐dependent protein deacetylase sirtuin‐1 SVA surrogate variable analysis SVM support vector machine SVM‐RFE support vector machine‐recursive feature elimination TFPI tissue factor pathway inhibitor TGFB1 transforming growth factor beta 1 TLR9 toll‐like receptor 9 VEGF vascular endothelial growth factor 1. Introduction Low back pain (LBP) has become the leading cause of disability worldwide, and it results in a loss of labor and enormous economic burden [[34]1]. Intervertebral disc degeneration (IVDD) is considered the main cause of LBP [[35]2]. However, conservative treatment is focused only on relieving symptoms; the histochemical changes that occur in IVDD cannot be reversed, and surgical therapy can cause other problems, including increased costs, adjacent segment degeneration, and trauma [[36]3, [37]4]. Thus, exploring the underlying mechanisms and identifying potential therapeutic strategies for IVDD are particularly important. The intervertebral disc (IVD), located between adjacent vertebrae of the spine, is composed of the central nucleus pulposus (NP), peripheral annulus fibrosus (AF) and cartilage endplate. The cells that compose the IVD include mainly NP cells, AF cells, and cartilage endplate cells. Among these cell types, NP cells can secrete extracellular matrix and relieve most of the mechanical pressure [[38]5]. It is believed that IVDD mainly occurs due to NP degeneration [[39]6]. Recently, the discovery of nucleus pulposus‐derived mesenchymal stem cells (NPMSCs) provides a novel direction for IVDD research. The endogenous repair ability of degenerative IVD may be mainly achieved by the differentiation of NPMSCs into NP cells [[40]7]. With increasing IVDD severity, the number of endogenous NPMSCs gradually decreases [[41]8, [42]9]. Wu et al. reported that NPMSCs derived from human degenerative IVD exhibit decreased cell proliferation and stem cell differentiation, decreased stem cell phenotypes, cell cycle arrest, and increased apoptosis; these findings indicate that the regenerative potential of NPMSCs decreases with the progression of IVDD [[43]10]. Our previous study also demonstrated that NPMSCs derived from degenerative IVD in rats display biological features of degeneration, such as reduced proliferation, impaired stemness maintenance, decreased expression of cell paracrine related genes hypoxia‐inducible factor‐1α (HIF‐1α), vascular endothelial growth factor (VEGF) and silent information regulator 1 (SIRT1) and increased cell apoptosis [[44]11]. Thus, IVDD may be caused by the decrease in the number and function of NPMSCs. In addition, the application of NPMSCs offers a novel strategy for treating IVDD. Mesenchymal stem cells (MSCs) have been extensively explored for the treatment of IVDD diseases [[45]12, [46]13, [47]14]. However, the microenvironment of IVDD is often unfavorable for the survival of exogenous MSCs [[48]15]. Compared with exogenous stem cells, NPMSCs are better adapted to the hypoxic and hyperosmotic microenvironments of degenerative IVD [[49]7, [50]11, [51]15]. Therefore, focusing research on NPMSCs is expected to reveal new mechanisms of IVDD and may provide perspectives for therapeutic development. Circadian rhythms usually refer to the internal system that temporarily regulates biological activity with an approximately 24‐h periodicity to allow an organism to cope with changes in the external environment [[52]16]. The circadian clock is located in the suprachiasmatic nucleus of the hypothalamus and transmits signals to peripheral tissues, synchronizing the peripheral biological clock. Circadian rhythms are controlled by a set of core genes, including brain and muscle aryl hydrocarbon receptor nuclear translocator‐like 1 (BMAL1), circadian locomotor output cycles kaput (CLOCK), period homologs 1, 2, and 3 (PER1/2/3), cryptochrome homologs 1 and 2 (CRY1/2), retinoid‐related orphan receptor (ROR), and nuclear receptor subfamily 1 group D (REV‐ERB). BMAL1 and CLOCK form heterodimers that bind to E‐BOX elements in the promoter regions of the PER and CRY genes to promote their transcription. When PER and CRY aggregate and reach a certain concentration in the cytoplasm, they form dimers and enter the nucleus to inhibit the transcriptional activity of BMAL1/CLOCK dimers, thereby inhibiting the transcription of PER and CRY themselves. In addition, the nuclear receptors ROR and REV‐ERB can bind to ROR‐binding elements in the promoter region of BMAL1 to promote and inhibit BMAL1 transcription, respectively. The nuclear receptor PPAR can bind to the PPRE‐binding element in the BMAL1 promoter region to promote BMAL1 transcription. These core clock genes regulate the rhythmic expression of hundreds of other genes and participate in complex physiological and pathological processes in the human body [[53]17]. Disruption of circadian rhythms can result in cancer, metabolic disorders, and neurodegenerative diseases [[54]18, [55]19, [56]20]. With increasing research, not only has the core rhythm gene BMAL1 been shown to be related to IVDD, but the clock control gene Nrf2 is also involved in IVDD [[57]21, [58]22]. However, the relationship between circadian rhythm‐related genes and IVDD remains unclear. The immune system also participates in IVDD, which occurs in three stages. First, IVDD damages the physical barrier, exposing the NP tissue to the immune system and inducing an initial inflammatory response, in which many immune cells and cytokines are involved [[59]23, [60]24]. Second, during self‐repair, granulation tissue gradually forms, and then vascular ingrowth exposes the NP tissue to immune cells further [[61]24]. Studies have shown that the presence of IgG antibodies targeting IVD‐specific proteins in human degenerative IVD and the extracellular matrix of NP tissues is considered an autoantigen that induces an autoimmune response [[62]23, [63]25]. Finally, due to the continuous immune cell activation and infiltration, released cytokines lead to an immune response similar to a positive feedback loop, accelerating the progression of IVDD [[64]26]. Circadian rhythm also regulates the immune system [[65]27]. For example, BMAL1:CLOCK heterodimer regulates toll‐like receptor 9 (TLR9) expression and represses the expression of the inflammatory monocyte chemokine ligand (CCL2) as well as REV‐ERBα suppressing the induction of interleukin‐6 [[66]28]. However, the relationship between immunity and rhythm in IVDD remains unclear. In the present study, we aimed to identify key differentially expressed rhythm‐related genes (DERGs) and further explore their relationship with the immune system through comprehensive bioinformatics analysis. Given the critical role of NPMSCs in IVDD, the expression and function of these DERGs were validated in NPMSCs. This study may be beneficial for exploring the potential molecular mechanisms underlying IVDD development and identifying possible molecular targets for IVDD treatment. 2. Methods and Materials 2.1. The Main Plan of the Study A flow diagram of this study was shown in Figure [67]1. Initially, two IVDD‐related gene expression profiles and rhythm‐related genes were obtained from the Gene Expression Omnibus and GeneCards databases. These genes were analyzed to identify DERGs. Analyses via Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and gene set enrichment analysis (GSEA) were conducted to explore the biological functions of these genes. Two machine learning algorithms were subsequently used to identify hub rhythm‐related genes. The subsequent steps involved analyzing immune cell infiltration and correlations with hub rhythm‐related genes, including immune factors, chemokines, immunoinhibitors, immunostimulators, MHC molecules, and receptors. Additionally, researchers have investigated the motif enrichment and predictive capabilities of miRNAs. Finally, NPMSCs were isolated to validate hub circadian rhythm‐related genes. FIGURE 1. FIGURE 1 [68]Open in a new tab Flow diagram presenting the main plan and process of the study. 2.2. Data Downloading The series matrix file for [69]GSE56081 [[70]29] was obtained from the NCBI GEO database, with the annotation file [71]GPL15314 covering the expression profiles of 10 patients‐5 controls and 5 with IVDD. Similarly, the series matrix file for [72]GSE70362 [[73]30] was sourced from the NCBI GEO database, with the annotation file [74]GPL17810, which documents expression data for 24 patients‐8 controls and 16 with IVDD (Table [75]1). The sequencing results of human NP tissue mRNAs were used for the following analysis. Batch effects were corrected via the SVA algorithm. Rhythm‐related genes with a relevance score above 1 were extracted from the GeneCards database ([76]https://www.genecards.org/) [[77]31]. TABLE 1. IVDD datasets included for analysis. Data set Platform information Control IVDD [78]GSE56081 [79]GPL15314 human NP tissue LncRNA and mRNA microarray 5 5 [80]GSE70362 [81]GPL17810 human NP and AF tissue mRNA microarray 8 16 [82]Open in a new tab 2.3. Screening Differentially Expressed Genes The Limma [[83]32] package, which is an R software tool for differential expression analysis, was used to identify DEGs among groups. Using the “Limma” R package, we analyzed the molecular mechanisms underlying IVDD and identified genes that were differentially expressed between the control and IVDD samples. The screening criteria for the DEGs were a p value < 0.05 and a |logFC| > 0.585. Specifically, volcano plots indicated the DEGs that were downregulated and upregulated. Heatmaps representing all the DEGs across all the samples were generated via the R package “heatmap”. 2.4. GO and KEGG Pathway Analysis The “ClusterProfiler” R package [[84]33] was used to annotate genes and explore the functional correlations among overlapping genes. The functional categories were assessed via GO [[85]34] and KEGG [[86]35] analyses. GO and KEGG pathway enrichment factors were considered significant if both the p value and the Q value were < 0.05. 2.5. Feature Selection Process of LASSO Regression and the SVM Algorithm The Lasso logistic regression and SVM [[87]36] algorithms were employed to identify biomarkers of the disease. The Lasso algorithm utilized the “glmnet” package. The optimal lambda (λ) value was determined using 10‐fold cross‐validation to minimize misclassification error and select the most relevant genes. Additionally, Support Vector Machine‐Recursive Feature Elimination (SVM‐RFE), a machine learning‐based feature selection method, was used to refine the selection of characteristic genes. This method iteratively removes the least significant features generated by the SVM model to optimize accuracy. The “e1071” software package was used to construct a model to assess the value of these biomarkers. 2.6. Enrichment Analysis of GSEA Pathways Predefined gene sets were subjected to GSEA to rank genes according to the degree of their differential expression between the control and IVDD samples, and then we checked whether the predefined gene sets were enriched at the top or bottom of the ranking table. GSEA facilitated the comparison of differences in KEGG analysis results between the high‐expression and low‐expression groups, enabling an investigation of the molecular mechanisms of key genes within these groups. The number of replacements was set to 1000, and the type of replacement was set to the phenotype. 2.7. Analysis of Immune Cell Infiltration The CIBERSORT method is extensively employed to assess various immune cell types within the microenvironment via gene signatures [[88]37]. Using support vector regression, we analyzed the expression matrix of immune cell subtypes through deconvolution. This method profiles 547 biomarkers, identifying 22 human immune cell phenotypes, such as T cells, B cells, plasma cells, and myeloid cell subsets. The CIBERSORT algorithm was used to analyze patient data, inferring the relative proportions of 22 types of infiltrating immune cells and examining the Pearson correlation between gene expression and immune cell counts. 2.8. Transcriptional Regulation Analysis of the Hub Genes The R package “RcisTarget” was used to predict transcription factors. RcisTarget is based on all the calculations of motifs. The normalized enrichment score (NES) for each motif was contingent upon the total number count of motifs in the database. In addition to the motifs that were annotated in the source data, we derived additional annotation files from motif similarity and gene sequence information. To estimate the overexpression of each motif within the gene set, we initially calculated the area under the curve (AUC) for each motif pair on the basis of their recovery curves in motif sequencing. The NES of each motif was calculated on the basis of the AUC distribution of all motifs in the gene set. 2.9. Isolation and Identification of NPMSCs Degenerative NP tissues were obtained from three patients with degenerative IVDD disease (Pfirrmann grades III‐V), whereas normal NP tissues were obtained from three patients with lumbar trauma (Pfirrmann grades I‐II). The study protocol was approved by the Ethics Committee of the Northern Jiangsu People's Hospital (Approval No. 202103427). Patient information was included in Table [89]S1 (Table [90]S1). The NP tissues were dissected under a dissecting microscope and then digested with 0.2% type II collagenase (Gibco, United States, catalog No. 17101015) for 8 h at 37°C in a 5% CO[2] atmosphere. The retrieved cells and partially digested tissues were subsequently washed twice with phosphate‐buffered saline (PBS), followed by centrifugation at 1000 rpm (rpm) for 3 min. The cells were then cultured at 1 × 10^5 cells/cm^2 in complete MSC medium (Cyagen, United States, catalog no. RASMX‐90011) at 37°C in a 5% CO[2], 5% O[2] environment. Specifically, cells were exposed to high‐concentration serum (50%) for 2 h, followed by incubation in standard culture medium [[91]38]. The culture medium was refreshed every 2 days, and the NPMSCs were maintained until they reached passage 3 before further experiments were performed. The cells were passaged at a ratio of 2:1, and the culture times varied depending on the experiment, but typically ranged from 2 to 3 days. Immunofluorescence staining was performed to assess the expression of surface markers that are associated with MSCs. The slides, which were 25 mm in diameter and coated with polylysine, were placed in a 12‐well plate, where the NPMSCs were seeded and cultured in complete MSC medium. Then, the cells were fixed with 4% paraformaldehyde for 15 min and rinsed twice with PBS supplemented with 0.5% Triton X‐100 for 15 min. Next, the cells were blocked with 10% bovine serum albumin for 1 h at 37°C and then incubated overnight at 4°C with primary antibodies against CD105 (Proteintech, China, catalog no. 10862‐1‐AP), CD90 (ABclonal, China, catalog no. A2126), CD73 (ABclonal, China, catalog no. A2029), CD34 (ABclonal, China, catalog no. A7429), Tie2 (Proteintech, China, catalog no. 19157‐1‐AP), and CD45 (ABclonal, China, catalog no. A2115) (1:100). The slides were subsequently washed twice with PBS and then incubated with secondary antibodies (Abcam, United Kingdom, catalog no. ab150077, ab150078) (1:500) for 1 h at room temperature. After treatment with antifade mounting medium containing 4',6‐diamidino‐2'‐phenylindole for 10 min, the cell slides were observed and images were captured via a fluorescence microscope (Leica, Wetzlar, Germany). Osteogenic, adipogenic, and chondrogenic differentiation were induced. NPMSCs were seeded in 6‐well plates and cultured until reaching approximately 80% confluency. Then, the culture medium was changed to osteogenic, cartilage, and adipogenic differentiation medium (OriCell, China, catalog No. RAXMX‐90041, RAXMX‐90021, RAXMX‐90031), and the medium was changed every 2 days. After reaching the deadline of induction, the cells were washed with PBS, and then the cells were processed with oil red O, alizarin red, and alcian blue. 2.10. Quantitative Real‐Time Polymerase Chain Reaction Total RNA was extracted with TRIzol reagent (Invitrogen, United States, catalog no. 15596–026). The reverse transcription of total RNA to synthesize complementary DNA (cDNA) and subsequent cDNA amplification were carried out following the manufacturer's instructions via a Prime Script‐RT reagent kit (Vazyme Biotech, China, catalog no. R123‐01) and SYBR Premix Ex Taq (Vazyme Biotech, China, catalog no. Q111‐02). Target gene expression in NPMSCs from different groups at different times was determined via the comparative Ct method. The data were normalized to the average gene expression of cells from the control group, which was used as the baseline for comparison. Specifically, the ΔΔCt method was applied, where GAPDH was used as the reference gene for internal normalization, and the control group was used for external normalization. The primers were designed on the basis of GenBank sequences via Prime 5.0 software and are detailed in Table [92]2. TABLE 2. Sequences of the primers used for real‐time PCR. Gene Primer sequence‐forward Primer sequence‐ reverse GAPDH 5′‐GGAGCGAGATCCCTCCAAAAT‐3′ 5′‐GGCTGTTGTCATACTTCTCATGG‐3′ CCND1 5′‐GCTGCGAAGTGGAAACCATC‐3′ 5′‐CCTCCTTCTGCACACATTTGAA‐3′ FOXO1 5′‐TCGTCATAATCTGTCCCTACACA‐3′ 5′‐CGGCTTCGGCTCTTAGCAAA‐3′ FRMD8 5′‐CGGGGCCAGAGTCTCTTTG‐3′ 5′‐TCGTGGCACTTGAGGAGATAG‐3 NTRK2 5′‐TCGTGGCATTTCCGAGATTGG‐3′ 5′‐TCGTCAGTTTGTTTCGGGTAAA‐3′ PRRT1 5′‐TCATCCGAAAAGTCAGGACTCC‐3′ 5′‐GTGCCAGACTGATGGTAGTGG‐3′ TFPI 5′‐TTGTGCATTCAAGGCGGATGA‐3′ 5′‐TCTTCGCACTGTCGAGTGAAA‐3′ HIF‐1α 5′‐GAACGTCGAAAAGAAAAGTCTCG‐3′ 5′‐CCTTATCAAGATGCGAACTCACA‐3′ SIRT1 5′‐TAGCCTTGTCAGATAAGGAAGGA‐3′ 5′‐ACAGCTTCACAGTCAACTTTGT‐3′ VEGF 5′‐AGGGCAGAATCATCACGAAGT‐3′ 5′‐AGGGTCTCGATTGGATGGCA‐3′ [93]Open in a new tab 2.11. Statistical Analysis All the statistical analyses were carried out in R (version 4.2.2) with Graph Pad Prism. All the statistical tests were bilateral, and p < 0.05 was considered to indicate statistical significance. The boxplot was generated using the R package ggplot2, and statistical significance was calculated using the Wilcoxon test. 3. Results 3.1. Identification, GO Analysis and KEGG Pathway Analysis of DERGs The [94]GSE56081 and [95]GSE70362 datasets related to IVDD were retrieved from the GEO database, and the data included expression profile data from 34 patients, including individuals in the control group (n = 13) and individuals in the IVDD group (n = 21). A PCA diagram was used to show the batch situation before and after correction. The findings indicated a reduction in the interchip batch effect (Figure [96]2a,b) following correction with the SVA algorithm. The limma package was used to compute the DEGs between the control and IVDD groups. A total of 372 DEGs (Figure [97]2c,d) were screened, including 240 upregulated genes and 132 downregulated genes. Then, rhythm‐related genes with a relevance score > 1 were extracted from the GeneCards database ([98]https://www.genecards.org/), and 44 overlapping genes (Figure [99]2e) were obtained. Pathway enrichment analysis of the overlapping genes revealed enrichment primarily in negative regulation of immune system processes, positive regulation of kinase activity, epithelial cell proliferation, and other pathways (Figure [100]3a). Additionally, KEGG enrichment analysis indicated enrichment primarily in the Hippo signaling pathway and the TGF‐beta signaling pathway (Figure [101]3b). FIGURE 2. FIGURE 2 [102]Open in a new tab Identification of DERGs. (a): PCA of gene expression in the GEO dataset before batch effect removal. (b): PCA of gene expression removal in the GEO dataset after batch effect. (c): Volcano plots of DEGs; x‐axis: Log2 (fold change); y‐axis: −log10 (adjusted p value); pink, black and blue nodes indicate that the DEGs that are upregulated, nonsignificantly different, and downregulated, respectively. (d): Heatmaps of DEGs; red: Degeneration group; blue: Control group. (e): Intersection map of the DEGs and rhythm‐related genes. FIGURE 3. FIGURE 3 [103]Open in a new tab GO and KEGG enrichment analyses of DERGs. (a): Bar chart of the GO enrichment analysis results, including biological process, cell component, and molecular function. (b): Bar chart of the results of the KEGG enrichment analysis. 3.2. Screening of Hub DERGs To further identify the key genes affecting IVDD, LASSO regression and SVM algorithms were used to screen the characteristic genes associated with IVDD. The LASSO regression model identified seven characteristic genes significantly associated with IVDD. The optimal regularization parameter (λ) was determined using 10‐fold cross‐validation, which minimized misclassification error (Figure [104]4a,b). The SVM‐RFE algorithm ranked genes based on their classification accuracy, selecting the top nine genes with the highest 5‐fold cross‐validation accuracy (94.3%) (Figure [105]4c). By intersecting the results from the two methods, six overlapping genes (Cyclin D1 (CCND1), Forkhead box O1 (FOXO1), FERM domain containing 8 (FRMD8), Neurotrophic tropomyosin receptor kinase b (NTRK2), Proline‐rich transmembrane protein 1 (PRRT1) and Tissue factor pathway inhibitor (TFPI)) were identified as hub genes (Figure [106]4d). These six genes, CCND1, FOXO1, FRMD8, NTRK2, PRRT1, and TFPI, were identified as key genes. The expression levels of the six hub genes were compared between normal and degenerative NP tissues. CCND1, FOXO1, FRMD8, NTRK2, and TFPI were significantly upregulated, whereas PRRT1 was significantly decreased in the IVDD group (Figure [107]4e). Moreover, these 6 key genes were inversely predicted with the miRcode database, and 86 miRNAs were found to be associated with the six key genes, and Cytoscape [[108]39] was used for visualization (Figure [109]4f). A list of these miRNAs was included in Table [110]S2. FIGURE 4. FIGURE 4 [111]Open in a new tab Identification of hub DERGs and the predicted mRNA–miRNA relationship pairs. (a)–(b): Hub DERGs was obtained via LASSO logistic regression model. (a) Coefficient profiles of genes across varying λ values. (b) Misclassification error rates for different λ values determined by 10‐fold cross‐validation, with the optimal λ highlighted. (c): Characteristic genes identified via the SVM‐RFE algorithm, showing the highest cross‐validation accuracy (0.943) for nine selected features. (d): Venn diagram illustrating the overlap between genes identified by the LASSO regression and SVM‐RFE algorithms. The overlap includes six hub genes (CCND1, FOXO1, FRMD8, NTRK2, PRRT1, and TFPI), indicating their consistency across different selection methods. (e): Boxplot showing the expression levels of six hub genes (CCND1, FOXO1, FRMD8, NTRK2, PRRT1, and TFPI) in normal (blue) and IVDD (pink) group. The values represent raw expression data. (f): The mRNA–miRNA relationship pairs of the hub DERGs. 3.3. Analysis of the Expression Levels of Core Clock Genes and Their Correlation With Hub Genes The hub genes were significantly correlated with multiple core clock genes (Figure [112]5a), among which FOXO1 was positively correlated with PER3 (Pearson r = 0.52) and PRRT1 was negatively correlated with PER3 (Pearson r = 0.602). The rhythm‐related genes were obtained from the GeneCards database ([113]https://www.genecards.org/). The differences in the expression levels of the first 12 genes in the Relevances core were analyzed, and the results revealed that there was no significant difference between the two groups (Figure [114]5b). FIGURE 5. FIGURE 5 [115]Open in a new tab Core clock genes in the control and IVDD samples. (a): Bubble plot of the relationships between the hub DERGs and core clock genes. (b): Box plot of the expression levels of core rhythm genes in the control and IVDD samples. 3.4. Enrichment Analysis of GSEA Pathways The specific signaling pathways associated with the hub genes were subsequently studied, and the potential molecular mechanism of the key genes affecting the progression of IVDD was explored. GSEA revealed that the pathways associated with high CCND1 expression were the pentose phosphate pathway and the ecm receptor interaction (Figure [116]6a). Moreover, the pathways associated with high FOXO1 expression were involved in sphingolipid metabolism, antigen processing, and presentation (Figure [117]6b). The pathways associated with high FRMD8 expression were the MAPK signaling pathway and B cell receptor signaling pathway (Figure [118]6c). The pathways associated with high NTRK2 expression were involved in oxidative phosphorylation and the cell cycle (Figure [119]6d). The pathways associated with high PRRT1 expression were involved in the T cell receptor signaling pathway and linoleic acid metabolism (Figure [120]6e). The pathways associated with high TFPI expression were the apoptosis and toll‐like receptor signaling pathways (Figure [121]6f). FIGURE 6. FIGURE 6 [122]Open in a new tab Line graph of GSEA results for the control and IVDD patient samples. (a): Pathways associated with high CCND1 expression. (b): Pathways associated with high FOXO1 expression. (c): Pathways associated with high FRMD8 expression. (d): Pathways associated with high NTRK2 expression. (e): Pathways associated with high expression of PPRT1. (f): Pathways associated with high TFPI expression. 3.5. Immune Cell Infiltration Analysis By analyzing the relationships between DERGs and immune cell infiltration in IVDD datasets, we further explored the potential molecular mechanism through which key rhythm‐related genes affect the progression of IVDD. The proportions of immune cells in each patient and the correlations between immune cells were shown in Figure [123]7a,b, which indicated that the top five by percentage are T cells CD4 memory resting, B cells memory, macrophages M0, macrophages M1, and Mast cells activated. The number of M0 macrophages in the IVDD group was significantly higher than that of the control group (Figure [124]7c). We further explored the relationships between key genes and immune cells, and the results revealed that several key genes were strongly correlated with immune cells (Figure [125]7d). Among them, M0 macrophages were significantly related to FRMD8, PRRT, and TFPI. T follicular helper cells were significantly related to FRDM8, FOXO1, and CCND1. FIGURE 7. FIGURE 7 [126]Open in a new tab Immune infiltration analysis. (a): Immune cell concentrations in each sample; the colors indicate immune cells. x‐axis: Id of patients. (b): Correlation analysis among immune cells; The heatmap shows the Pearson correlation coefficients between different immune cell types. The colors of the cells represent the strength of the correlation, with red indicating a strong positive correlation and blue indicating a strong negative correlation. (c): Violin plot of the difference in each type of immune cell between the control and IVDD samples; x‐axis: Immune cells; y‐axis: Cell concentration; red: IVDD samples; blue: Control samples. (d): Correlation analysis between the hub DERGs and various infiltrating immune cells; x‐axis: Immune cells; y‐axis: Hub DERGs. The colors of the cells represent the strength of the correlation, with red indicating a positive correlation and blue indicating a negative correlation. The color intensity of each cell in the heatmap corresponds to the magnitude of the correlation, with darker shades indicating stronger correlations. 3.6. Correlation Analysis Between Hub DERGs and Immune Factors The correlations between the 6 hub genes and different immune factors, including immunosuppressive factors, immunostimulatory factors, chemokines, MHCs, and receptors, were analyzed. CCND1 was positively correlated with TAP1 (Figure [127]8d). FOXO1 was positively correlated with HLA‐DRA, CCR7, and XCR1 but negatively correlated with CCL13 (Figure [128]8a,d,e). RND8 was positively correlated with TGFB1, LTA, HLA‐DRA, and XCR1 but negatively correlated with CCL13, CXCL6, and TNFRSF25 (Figure [129]8a–e). NTRK2 was positively correlated with TGFB1 and IL10RB but negatively correlated with CCL13, CXCL6, and KDR (Figure [130]8a,b). PRRT1 was positively correlated with CXCL12, TNFRSF8, CD48, and XCR1 but negatively correlated with LTA, HLA‐DRA, and CCR7 (Figure [131]8a–e). TFPI was positively correlated with CXCL8 and TGFB1 (Figure [132]8a,b). These analyses revealed that key genes are closely related to the level of immune cell infiltration and play important roles in the immune microenvironment (Figure [133]8a–e). FIGURE 8. FIGURE 8 [134]Open in a new tab Correlation analysis between hub DERGs and immune factors. (a): Correlation analysis between hub DERGs and chemokines. (b): Correlation analysis between hub DERGs and immunoinhibitors. (c): Correlation analysis between hub DERGs and immunostimulators. (d): Correlation analysis between hub DERGs and MHCs. (e): Correlation analysis between hub DERGs and receptors. 3.7. Analysis of the Transcriptional Regulation of Key Genes The 6 hub DERGs (CCND1, FOXO1, FRMD8, NTRK2, PRRT1, and TFPI) described above in this analysis were used, and the results revealed that they were regulated by numerous transcription factors and other common mechanisms. Estimates of the global mean and standard deviation for the hub DERGs were shown in Figure [135]9a. Therefore, these transcription factors were enriched and further analyzed via a cumulative recovery curve. The transcription factors that regulate these genes were identified using the “motifAnnotations_hgnc_v9” dataset, and their binding motifs were further analyzed using a cumulative recovery curve. The results of motif‐TF annotation and selection analysis of important genes revealed that the motif with the highest standardized enrichment score (NES: 5.37) was cisbp_M2567, which was enriched in the FOXO1, FRMD8, NTRK2, and PRRT1 (Figure [136]9b,e). The motif with the second standardized enrichment score (NES: 4.46) was cisbp_M0619, which was enriched in the FOXO1, FRMD8, and NTRK2 (Figure [137]9c,e). The motif with the third standardized enrichment score (NES: 4.31) was cisbp_M0413, which was enriched in the FOXO1, FRMD8, and NTRK2 (Figure [138]9d,e). The binding motifs were shown in the first line of Figure [139]9e. FIGURE 9. FIGURE 9 [140]Open in a new tab Transcriptional regulation analysis of CCND1, FOXO1, FRMD8, NTRK2, PRRT1 and TFPI. (a): Estimates of the global mean and standard deviation for all the hub DERGs. (b–d) The enrichment analysis results of different transcription factor binding sites associated with the hub DERGs. (e): Major transcription factors and binding sites of the hub DERGs, with the normalized enrichment scores (NES) shown for each motif. 3.8. Identification of NPMSCs and the Expression of Hub DERGs in NPMSCs at Different Time The MRI images of the IVD tissue were shown in Figure [141]10a. Cells isolated from the normal IVD tissue were used as the control group, and the cells from the degeneration IVD tissue were used as the degeneration group. The cells isolated from the human IVD presented a long spindle shape. As shown in Figure [142]10b, the cells exhibited strong positive staining for MSCs markers CD90, CD105, and CD73, as well as the NP cell marker Tie2, confirming their mesenchymal and NP origin. NPMSCs were positive for Alizarin red, Oil Red O, and Alcian blue staining after induced differentiation, and NPMSCs derived from degenerative tissues exhibited reduced differentiation potential compared to the control group (Figure [143]10c). FIGURE 10. FIGURE 10 [144]Open in a new tab Identification NPMSCs and the expression of hub DERGs in NPMSCs. (a): MRI results of the samples. (b): Immunofluorescence analysis of MSC‐associated surface markers (250 μm). (c):Trilineage differentiation of NPMSCs from degeneration and normal tissues. (d): Expression of HIF‐1α, VEGF and SIRT1 in normal and degenerated NPMSCs. (e): CCND1, FOXO1, FRMD8, NTRK2, PRRT1 and TFPI expression in normal and degenerated NPMSCs at different time. The data are presented as the mean ± SD, n = 3; *p < 0.05 vs. the control group; as evaluated by the t test, ZT, Zeitgebers (German for “time giver”). Compared with the control group, the degeneration indices, including HIF‐1α, VEGF, and SIRT1, were decreased in the degeneration group (Figure [145]10d). This result indicated that NPMSCs derived from degenerative IVD exhibit degenerative characteristics [[146]9, [147]11, [148]40, [149]41]. The expression levels of FOXO1 followed a 12‐h rhythm in both the control and degeneration groups. Throughout the rhythm, FOXO1 showed a higher level of expression compared to the control group. The expression levels of CCND1, FRMD8, NTRK2, and TFPI were dysrhythmic in the degeneration group compared to the control group, exhibiting a 12‐h rhythm in some cases. However, PRRT1 did not seem to show a significant circadian rhythm in either group (Figure [150]10e). 4. Discussion As the main cause of LBP, IVDD is a common degenerative disease that leads to physical pain and economic burdens on societies worldwide. However, the underlying pathological mechanism is still unclear, and currently available treatment measures cannot reverse the pathological changes that occur in IVDD. Therefore, identifying potential targets and possible effective drugs for the treatment of IVDD will lead to the development of new approaches for the prevention and treatment of LBP. With the development of society, circadian rhythm disorders and related diseases caused by cross‐time flight, frequent shifts, and exposure to light at nighttime are becoming increasingly common. Studies have shown that circadian rhythms are associated with IVDD [[151]42], and several factors, such as oxidative stress and mechanical pressure, can also lead to rhythm disorders in IVDD. Excessive loading can disrupt the circadian rhythm of the NP via the RhoA/ROCK pathway, which leads to IVDD accompanied by disrupted rhythmic gene expression, and reversing the expression of BMAL1 can protect against compression‐induced IVDD [[152]21]. In BMAL1 knockout mice, the loss of circadian rhythms resulted in age‐related changes in extracellular matrix (ECM) structure and an endochondral ossification‐like phenotype in AF cells, with dysregulation of ECM and FOXO pathways. This was further associated with FOXO1 upregulation, which impaired osteogenic differentiation in AF cells, suggesting potential molecular targets for IVDD [[153]43]. Interestingly, recent research found that MSCs, which have potential for regenerative therapies in IVDD, also exhibit circadian rhythm regulation of core rhythm genes [[154]15, [155]44]. BMAL1 deletion in joint MSCs led to disrupted circadian rhythms, increased pro‐inflammatory responses, and exacerbated arthritis, emphasizing the key role of core rhythm genes in maintaining joint health and regulating inflammatory responses [[156]44]. However, core rhythm genes did not show significant differences in the dataset in this study. Nevertheless, it remains meaningful to conduct further experimental validation in NPMSCs in future studies. In addition to the core rhythm gene, the clock‐controlled gene nuclear factor erythroid 2‐related factor 2 (NRF2) can also alleviate IVDD by alleviating inflammation‐induced rhythm disorder of NP cells [[157]22]. These findings revealed that there is a network of mechanisms that link rhythm genes and IVDD. Thus, several DERGs can be used as molecular targets for treatment. The development of bioinformatics has accelerated investigations on the internal mechanism of disease. Through bioinformatics analysis, the six hub DERGs were obtained. Then, through experimental verification, we found that the expression levels of CCND1, FRMD8, and NTRK2 exhibited a more stable pattern in the degeneration group compared to the control group, suggesting that their circadian rhythmicity may be suppressed or altered during the degeneration process. In contrast, TFPI showed no significant circadian oscillation in the degeneration group, indicating that its rhythmicity may be disrupted. Additionally, the FOXO1 expression level was upregulated in the degeneration group. These results indicated that some DERGs in degenerate NPMSCs lose control of circadian rhythm compared to normal NPMSCs, but some still maintain rhythmicity, which needs further study. According to several bioinformatics studies, CCND1 and NTRK2 may be of diagnostic importance and are related to IVDD [[158]45, [159]46]. CCND1 and FOXO1 have been proven to play essential roles in IVDD [[160]47, [161]48, [162]49], whereas FRMD8, PRRT1, and TFPI have not been confirmed to be associated with IVDD. In addition, functional enrichment analysis indicated that these hub genes were involved in immune function, metabolism, cell cycle, and death. Hence, the results of the present study may provide new insight into the diagnosis and treatment of IVDD. The number and function of NPMSCs, which are an important component of NP tissue, are decreased in degenerative IVD [[163]50]. Our previous study revealed that NPMSCs derived from degenerative IVD exhibit evidence of degeneration [[164]11, [165]41]. However, the underlying mechanism is unclear. Thus, exploring this process was particularly important. In the present study, the expression of DERGs was verified in NPMSCs, and these cells were used to reveal the relationship between rhythm genes and NPMSCs degeneration. The results revealed that DERGs are differentially expressed between degenerative and normal NPMSCs, indicating that rhythm is related to NPMSCs degeneration. However, the specific mechanism still needs to be further studied. CCND1, which is a D‐type cyclin protein (Dcyclin), is a key element in the control of cell cycle progression [[166]51] and contributes to the proliferation of NP cells. Cyclin D1 is significantly upregulated in a mechanically induced IVDD model [[167]47, [168]49]. The decrease in CLOCK abundance in HC11 cells was related to the increase in CCND1 expression, which was consistent with the increase in the growth rate and cell cycle progression of shCLOCK‐transfected cells [[169]52]. Thus, CCND1 may have an effect on IVDD by modulating rhythm. However, the expression of CCND1 in megakaryocytes does not seem to exhibit obvious rhythmicity [[170]53], and this study revealed that CCND1 is not closely related to the core rhythm genes. Therefore, the rhythmicity of CCND1 in IVD needs to be verified by further study. The pathways associated with high CCND1 expression were the pentose phosphate pathway and ECM–receptor interaction. The pentose phosphate pathway plays a crucial role in cellular metabolism and oxidative stress responses, which are critical in IVDD [[171]54]. ECM–receptor interactions are also central to disc degeneration, where changes in extracellular matrix remodeling contribute to the pathological process of IVDD [[172]55]. FOXO1 protein is a transcription factor that performs important functions in development, senescence, and longevity, and its expression is significantly downregulated in degenerative IVD [[173]56]. Together with SIRT1, FOXO1 participates in IVDD through oxidative stress and has certain therapeutic significance [[174]57]. In addition, FOXO1 may be activated during deficiency of the PH domain leucine‐ richrepeat protein phosphatase (PHLPP1) to protect against NP cell phenotypic changes, extracellular matrix degradation, and cell apoptosis during IVDD [[175]58]. A deficiency in FOXO1 can impair the homeostasis of the IVD, affect the differentiation and maturation of NP cells, and regulate autophagy and antioxidative processes in NP cells [[176]48]. An imbalance in the FOXO1 pathway in patients with osteoarthritis can lead to changes in the circadian rhythm and dynamic balance of the extracellular matrix [[177]59]. These studies revealed that FOXO1 may also be rhythmically regulated in IVDD, and this was confirmed in the present study. In the present study, FOXO1 was linked to sphingolipid metabolism and antigen processing and presentation. Sphingolipid metabolism is involved in the inflammatory response [[178]60], and antigen processing and presentation suggest a role for immune modulation [[179]61], which supports our hypothesis that immune dysfunction is an important factor in IVDD progression. In addition, FOXO1 was found to be closely related to the core clock genes PER1, PER3, CRY1, and follicular helper T cells in IVDD. Therefore, together with these core clock genes, FOXO1 may participate in IVDD through the senescence pathway and immune system. NTRK2 is a member of a family of transmembrane tyrosine kinases that were previously thought to be related to neurogenesis, but it has also been shown to be related to the occurrence and development of many diseases, such as cancer, kidney disease, and osteoarthritis [[180]62, [181]63, [182]64]. Pathways linked to NTRK2 expression include oxidative phosphorylation and cell cycle regulation. These pathways are crucial for energy metabolism and cell division, which could contribute to tissue remodeling in degenerating intervertebral discs [[183]65, [184]66]. There is a significantly positive correlation between the mRNA expression of NTRK2 and the expression of the immune chemokine CX3CL1 in osteoarthritic synovia [[185]62]. However, according to the results of the present study, there was no significant correlation between these two parameters in IVDD. CRY1 was found to be closely related to NTRK2 in the present study, which showed that NTRK2, along with CRY2, may participate in IVDD. TFPI is a multivalent Kunitz‐type serine protease inhibitor that is involved in extracellular matrix degradation and remodeling in tumors [[186]67]. IVDD is also characterized by extracellular matrix degradation and remodeling, and an investigation of temporal fluctuations in constant darkness indicated the existence of daily rhythms for TFPI [[187]68], which suggests that the circadian rhythm may be involved in IVDD through TFPI. The pathways associated with TFPI were apoptosis and toll‐like receptor signaling pathways. Both apoptosis and toll‐like receptor pathways are involved in the inflammatory processes and cell death that contribute to IVDD progression [[188]69]. We found that TFPI was significantly correlated with the numbers of resting CD4^+ memory T cells and M0 macrophages. TFPI can affect macrophage apoptosis and alleviate inflammation induced by endoplasmic reticulum stress in human M2‐polarized macrophages [[189]70, [190]71], suggesting that a similar mechanism may exist in IVDD, which needs to be further studied later. PRRT1 is a component of native AMPAR complexes in multiple regions of the brain, and its deletion can lead to altered surface levels and phosphorylation of AMPARs as well as impaired forms of synaptic plasticity [[191]72]. The present study found that PRRT1 is positively correlated with PER1 and CRY1, negatively correlated with CLOCK, DBP, and PER3, and closely related to plasma cells. The pathways enriched in PRRT1 were the T cell receptor signaling pathway and linoleic acid metabolism. The T cell receptor signaling pathway pointed to immune cell activation, further supporting the role of immune responses in IVDD. Linoleic acid metabolism is also implicated in inflammation, which is a known factor in IVDD [[192]73]. FRMD8 can promote inflammation and growth factor signal transduction by stabilizing the iRhom/ADAM17 exfoliase complex in human macrophages [[193]74]. The present study revealed that PRRT1 is also related to macrophages and CLOCK, indicating that PRRT1 may interact with CLOCK and affect IVDD together with M0 macrophages. Although the roles of M1, M2, and Th macrophages in IVDD have been studied, the roles of M0 macrophages and follicular helper T cells, as well as their underlying mechanisms, in IVDD have not yet been revealed. Our analysis revealed that M0 macrophages and follicular helper T cells are most closely related to the six hub genes during IVDD. Previous studies have shown that macrophages are the only inflammatory cells that can penetrate the NP, and the number of macrophages is positively correlated with the severity of IVDD [[194]75]. Macrophages participate in all stages of IVDD, and the proportion of M2 macrophages decreases with the progression of IVDD, whereas the proportion of M1 macrophages increases [[195]76]. Yao et al. reported that the number of Th1 cells increased significantly on the 14th day but decreased on the 28th day, and the number of Th2 cells increased on the 28th day after the establishment of an IVDD rat model, which was synchronized with the transformation of macrophages from the M1 phenotype to the M2 phenotype in the IVD [[196]77]. Geiss et al. reported that NP exposure via the immune system for 3 weeks can promote the development of T cells into Th2 cells and exert their effective function by releasing IL‐4 [[197]78]. Thus, the concentrations of macrophages and Th cells vary with the progression of IVDD, and these cells participate mainly in IVDD through inflammatory factors. But the exact role and mechanism of M0 macrophages and follicular helper T cells in IVDD need to be further explored. The regulatory mechanisms of circadian rhythm genes in IVDD remain largely unexplored. In this study, we performed motif enrichment analysis to identify potential transcription factors (TFs) that may regulate the six hub genes (CCND1, FOXO1, FRMD8, NTRK2, PRRT1, TFPI) in IVDD. The cisbp_M2567 motif, which is enriched in FOXO1, FRMD8, NTRK2, and PRRT1, is recognized by heat shock transcription factors (HSF1, HSF2, HSF4), suggesting that these genes may be involved in the cellular stress response. HSF1 exhibits increased expression in herniated discs, suggesting a protective stress response to environmental factors associated with disc degeneration [[198]79]. Moreover, ZBTB1, a transcription factor associated with the cisbp_M0413 motif, was found to regulate FOXO1, FRMD8, and NTRK2, indicating a possible link between circadian gene dysregulation and immune dysfunction. ZBTB1 prevents DNA damage in replicating immune progenitors, allowing the generation of B cells, T cells, and myeloid cells [[199]80], further supporting the hypothesis that immune dysregulation contributes to IVDD pathology. Additionally, the OVOL1/OVOL2/OVOL3 transcription factors were identified as regulators of FOXO1 and PRRT1, highlighting their potential role in cell differentiation and metabolic homeostasis [[200]81]. These findings suggest that the rhythmic disruption of these hub genes in IVDD is not only driven by circadian clock dysregulation but also by transcriptional alterations mediated by stress response and immune‐related transcription factors. Future studies should investigate whether modulating these TFs can restore normal circadian rhythm and mitigate IVDD progression. The present study also has several limitations. First, only a few DNA microarray datasets of IVDD are available in the GEO database, and clinical validation in large samples is lacking. Second, the exact molecular mechanism and circadian rhythmicity of the hub genes involved in IVDD still need to be further verified via in vivo and in vitro experiments. Additionally, it is important to note that the observation period (24 h) may be insufficient to conclusively determine the presence or absence of a 24‐h rhythm in gene expression. Longer observation times, spanning multiple circadian cycles, are necessary to confirm the rhythmicity and to distinguish between 12 and 24‐h cycles. It would be of interest to investigate whether core clock genes are expressed in cyclic modes in NPMSCs, as understanding their rhythmic expression could provide deeper insights into the rhythm‐related molecular mechanisms in NPMSCs. Finally, the retrospective nature of the research limits the clinical value of this work, indicating that multicenter or prospective studies are imperative for elucidating the relationship between rhythm and IVDD. 5. Conclusion In conclusion, our research provides innovative perspectives on the role of rhythms in IVDD. CCND1, FOXO1, FRMD8, NTRK2, PRRT1, and TFPI were successfully demonstrated to be hub rhythm‐related genes that, along with immunity, are involved in the procession of IVDD. This study may provide new ideas for elucidating the molecular mechanism and identifying potential targets for treating IVDD. Author Contributions All the authors contributed to the study conception and design. Data curation, Yongbo Zhang, Sheng Yang, and Liuyang Chen; formal analysis, Yongbo Zhang and Sheng Yang; funding acquisition, Liang Zhang; methodology, Sheng Yang; project administration, Liang Zhang; validation, Liuyang Chen; visualization, Liang Zhang and Rui Dai; writing‐original draft, Yongbo Zhang, Sheng Yang, Liuyang Chen, and Hua Sun; writing‐review and editing, Yongbo Zhang, Liang Zhang, Rui Dai, and Sheng Yang. All the authors drafted and reviewed the manuscript. Ethics Statement The experimental protocol was established according to the ethical guidelines of the Helsinki Declaration and was approved by the Human Ethics Committee of Northern Jiangsu People's Hospital (Approval No. 202103427). Written informed consent was obtained from all participants. Clinical trial number: not applicable. Conflicts of Interest The authors declare no conflicts of interest. Supporting information Table S1. Patient information summary. Table S2. mRNA–miRNA relationship pairs. [201]JSP2-8-e70066-s001.docx^ (20.8KB, docx) Acknowledgments