Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Med., 12 November 2025

Sec. Translational Medicine

Volume 12 - 2025 | https://doi.org/10.3389/fmed.2025.1698050

This article is part of the Research TopicThe Application of Multi-omics Analysis in Translational MedicineView all 15 articles

Identification and validation of novel risk genes for intervertebral disc disorder by integrating large-scale multi-omics analyses and experimental studies


Zhe Zhang,Zhe Zhang1,2Yu ChenYu Chen3Qianjin WangQianjin Wang1Zheng LiZheng Li4Bingyang Dai,Bingyang Dai5,6Cheng Dong,Cheng Dong5,6Zhengya Zhu*Zhengya Zhu2*
  • 1Department of Orthopaedics and Traumatology, Faculty of Medicine, The Chinese University of Hong Kong, Hong Kong SAR, China
  • 2Department of Orthopaedics, The Affiliated Hospital of Xuzhou Medical University, Xuzhou, Jiangsu, China
  • 3Division of Spine Surgery, Department of Orthopaedic Surgery, Nanjing Drum Tower Hospital, Affiliated Hospital of Medical School, Nanjing University, Nanjing, China
  • 4Department of Orthopaedics, The Affiliated Hospital, Southwest Medical University, Luzhou, Sichuan, China
  • 5Department of Biomedical Engineering, Faculty of Engineering, The Hong Kong Polytechnic University, Hong Kong, Hong Kong SAR, China
  • 6The Hong Kong Polytechnic University Shenzhen Research Institute, Shenzhen, China

Introduction: Although genome-wide association studies (GWAS) have identified multiple genetic loci linked to intervertebral disc disorder (IDD), their functional characterization remains largely unelucidated. We aim to leverage an integrative analytical pipeline to identify novel IDD risk genes from genetic associations and experimentally validate their functional roles.

Methods: We integrated transcriptome-wide association studies (TWAS), proteome-wide association studies (PWAS), expression and protein quantitative trait loci (eQTL and pQTL) colocalization analyses to identify potential causal genes for IDD. Enrichment analysis, expression profiling, protein-protein interaction (PPI) network construction, and druggability evaluation were also performed for the prioritized causal candidates. Subsequently, human intervertebral disc (IVD) tissues spanning degeneration grades and an in vivo mouse IDD model were employed to functionally characterize candidate risk genes.

Results: Integrative analysis of TWAS and PWAS with colocalization studies identified 104 genes and 10 proteins exhibiting causal associations with IDD. The identified genes/proteins were enriched in extracellular matrix organization, cellular senescence and collagen formation. Crucially, TMEM190, CILP2, and FOXO3 were demonstrated consistent evidence across TWAS, two independent PWAS datasets, and corresponding colocalization analyses, with CILP2 emerging as a potentially druggable target. Differential expression analysis revealed significant upregulated TMEM190 and CILP2, along with downregulated FOXO3 during IVD degeneration. These results were subsequently confirmed at protein levels in clinical specimens. Mouse model experiments further established that down-regulation of CILP2 alleviated IDD progression.

Discussion: Collectively, this work provides an updated compendium of putative IDD risk genes and delineates pathogenic roles for TMEM190, CILP2, and FOXO3, providing a broad hint for further research on novel mechanism and therapeutic targets for IDD.

Introduction

Intervertebral disc disorder (IDD) represents a primary etiology of low back pain (1). Both genetic predispositions and environmental risk factors involved in its pathogenesis (25). Currently, the treatment of IDD primarily relies on symptomatic management with NSAIDs and surgical interventions for more severe cases. However, symptomatic treatments fail to address underlying disease mechanisms, while surgery entails significant costs, potential complications, and surgical morbidity (6). Therefore, identifying causative genes and developing targeted therapeutic strategies is imperative for IDD management.

Recent GWAS have identified multiple loci associated with IDD, predominantly within non-coding genomic regions (7, 8). These regions exhibit complex regulatory mechanisms and linkage disequilibrium, complicating the identification of underlying causal genes. TWAS coupled with eQTL colocalization address this limitation by linking non-coding disease-associated variants to transcriptional changes. In a TWAS study, genetic predictors of gene expression, specifically cis-eQTLs regulating nearby genes, are identified in reference populations, such as the Genotype-Tissue Expression (GTEx) project. These genetic predictors subsequently impute transcriptomic profiles in GWAS cohorts to identify gene expression levels and disease traits (9). eQTL colocalization analyses determine whether shared causal variants gene expression and disease risk share the same causal variants underlie both gene expression and disease risk, strengthening causal inference for candidate genes (10). While IDD-specific TWAS study remains scarce due to the limited large-scale human transcriptomic datasets of disc tissues, GTEx demonstrates substantial eQTLs conservation across tissues (11, 12). Thus, regulatory variants identified in non-disc tissues may still modulate disc biology and IDD susceptibility.

Complementary to TWAS, PWAS utilizes pQTL data to identify protein-level associations with diseases, providing enhanced mechanistic insight. Recently, large-scale human plasma proteome datasets, including those from the ARIC study and Iceland Biobank, have enabled robust pQTL derivation, facilitating practical PWAS implementation (13). Plasma proteins, which serve as key druggable targets and biomarkers for complex traits, can reflect systemic pathological changes associated with IDD. While PWAS has been applied to other diseases (1416), its application to IDD remains unexplored. Future integration of PWAS with TWAS and eQTL/pQTL colocalization will enable the identification of disease-causing genes with higher precision, while minimizing confounding effects from horizontal pleiotropy. This multimodal approach will elucidate IDD molecular mechanisms and accelerate development of targeted therapeutics.

This study aimed to identify and validate potential causal genes associated with IDD by integrating multi-omics analyses with experimental approaches. We performed both TWAS and PWAS to uncover novel genes implicated in the pathogenesis of IDD (Figure 1). Colocalization analyses were performed to establish potential causal relationships between these genes and IDD risk. The expression patterns of the prioritized genes in human degenerated intervertebral discs were also assessed. Furthermore, enrichment analysis was conducted to identify key pathways and biological terms associated with IDD. Additionally, we explored protein-protein interactions among the candidates and evaluated their druggability. Finally, experimental validation of top-prioritized genes (TMEM190, CILP2, FOXO3) was conducted using clinical specimens and animal model.

FIGURE 1
Flowchart depicting the integration of transcription prediction models, protein level prediction models, and GWAS summary statistics for intervertebral disc disorder (IDD). It begins with predicting gene expression and proteins, identifying associated genes and proteins for IDD, and performing colocalization analysis. Cross-validation identifies three candidate genes: TMEM190, CILP2, and FOXO3. Validation includes clinical and animal studies, enrichment, gene expression, protein interaction, and drug analysis.

Figure 1. Overview of the study. This schematic illustrates the multistep approach employed to identify potential genes associated with IDD. First, TWAS and eQTL analyses were conducted to identify potential risk genes. Second, two independent PWAS and pQTL analyses were performed to identify potential causal proteins of IDD. Third, enrichment analyses were conducted to elucidate the functions of these potential causal genes/proteins. Fourth, data from the GEO were analyzed to identify differentially expressed potential causal genes. Additionally, PPI analyses were performed to explore interactions among the identified genes. Furthermore, the druggability of the potential causal genes and available drugs that target these genes were explored. Finally, validation studies were performed with clinical samples and an animal model.

Materials and methods

IDD GWAS summary data sources

The GWAS data were obtained from meta-analyses of data from several large cohorts, including deCODE Genetics from Iceland, the Danish Blood Donor Study, the Copenhagen Hospital Biobank, and the UK Biobank (7). Participants provided blood or buccal samples with informed consent, permitting the use of their samples and data in deCODE Genetics and the UK Biobank dataset were included. Each dataset was assumed to share a common odds ratio (OR), while allowing for different population frequencies of alleles and genotypes. Variants with imputation information scores below 0.8 were excluded from the analyses. The GWAS summary data used in our analysis come from the worldcome y data used in These analyses included 58,854 IDD cases and 922,958 controls, with the participants being of European descent. A total of 53.5 million sequence variants were included in the GWAS analysis.

Transcriptomic data from multiple human tissues

The eQTL data were obtained from GTEx Version 8 (49 tissues) (11). GTEx provides extensive data on the relationship between genetic variation and gene expression, sourced from 838 postmortem donors, and 15,201 RNA-sequencing samples were included, primarily of European descent.

Human protein abundance references in discovery proteome-wide association studies

The pQTL datasets incorporated in this study were derived from two large-scale investigations: the ARIC study, which includes data on 4,423 proteins from 7,213 individuals (13), and deCODE Genetics, which encompasses 4,428 proteins from 35,559 individuals (17), both primarily of European descent.

Human intervertebral disc degeneration microarray datasets

Human disc tissue expression microarray datasets were obtained retrospectively from the GEO (GSE56081: n = 10, with five degenerative and five normal samples).

Transcriptome-wide association studies

We performed TWAS analysis by integrating genome-wide summary statistics from an IDD GWAS with eQTL data from GTEx Version 8 across 49 tissue types as descripted before (18). To ensure consistency between datasets, we harmonized the IDD GWAS single nucleotide polymorphisms (SNPs) with the GTEx reference data, aligning SNP reference alleles, effect alleles, and associated metadata. Single-tissue TWAS was conducted for all tissues via SPrediXcan, followed by cross-tissue analysis via S-MultiXcan. S-MultiX can, which is based on a multitissue integration approach, allows for the combination of gene expression data across tissues, enhancing statistical power and enabling the identification of candidate susceptibility genes. We utilized default parameters for the software, with the exception of adjusting the “—cutoff_condition_number” parameter to 30. Only protein-coding genes were considered in the analysis, and significance was determined via a false discovery rate (FDR) threshold of p < 0.05.

FastENLOC colocalization

FastENLOC colocalization tool was used to strengthen the causal inferences drawn from our TWAS findings (19). Briefly, we computed posterior inclusion probability (PIP) values from IDD GWAS data via the torus tool, which quantifies the likelihood of each SNP’s association with IDD. These PIPs were then input into fastENLOC, which performs colocalization analysis by integrating GWAS PIP values with precomputed GTEx multitissue eQTL annotations. Colocalization was performed for each tissue, producing gene-level colocalization probabilities (GLCPs), indicating the likelihood that a given variant influences both IDD GWAS and gene expression in each tissue. The results across all the tissues were then merged, and for each gene, the maximum GLCP value was retained to identify the tissue with the strongest colocalization signal. Genes with a maximum GLCP > 0.5 were considered to have significant evidence of colocalization.

Proteome-wide association studies

BLISS method was used for PWAS analysis (20). Traditional PWAS approaches rely on individual-level proteomic data, which restricts the use of extensive summary-level pQTL datasets available in public repositories. The BLISS method enables the utilization of large-scale summary-level datasets for more efficient proteomic association analysis by constructing protein imputation models directly from summary-level pQTL data. In this study, we performed PWAS analyses via summary-level pQTL data from two large-scale cohorts: the ARIC study and deCODE Genetics. These datasets include over 4,000 proteins, facilitating a comprehensive analysis of protein–trait associations in the context of IDD. For discovery purposes, proteins with a nominal p < 0.05 were considered significant.

Bayesian colocalization analysis

We also performed Bayesian colocalization analyses via the coloc R package to investigate whether the identified associations between plasma proteins and IDD share the same causal variants rather than being affected by linkage disequilibrium (21). The Bayesian colocalization method evaluated evidence for five distinct hypotheses at each locus: (1) no association with either trait, (2) association with trait 1 only, (3) association with trait 2 only, (4) both traits are associated, but each has distinct causal variants, and (5) both traits share a common causal variant (22). Posterior probabilities for each hypothesis (H0, H1, H2, H3, and H4) were calculated as part of the analysis. Initial prior probabilities were assigned as follows: a SNP exclusively associated with trait 1 (p1) had a probability of 1 × 10−4, a SNP exclusively associated with trait 2 (p2) had a probability of 1 × 10−4, and a SNP shared by both traits (p12) had a probability of 1 × 10−5 (23). A posterior probability of H4 (PPH4) > 0.5 was considered evidence of a shared causal variant between the two traits.

Enrichment analysis of significant findings

To further investigate the biological role of the identified genes, we performed enrichment analysis via Gene Ontology (GO) categories [encompassing biological processes (BPs), molecular functions (MFs), and cellular components (CCs)], Kyoto Encyclopedia Genes and Genomes (KEGGs) pathways, and Reactome pathways (24, 25). The analysis was performed via the clusterProfiler and Reactome PA R packages. Significant genes or proteins were defined as those with a p < 0.05. The background set for the enrichment analysis consisted of all genes or proteins tested in the study, representing the total gene/protein pool from which the significant findings were derived. The ggplot2 R package was used for visualization.

Annotation of prioritized genes/proteins

The genes/proteins prioritized through TWAS, PWAS, and colocalization were further annotated by evaluating their expression levels in degenerative disc tissues and constructing gene coexpression networks. First, we obtained human disc tissue expression microarray datasets from the GEO (GSE56081: n = 10, with five degenerative and five normal samples).1 After normalization of the expression matrix, differential expression analysis was performed via the lmFit() and eBayes() functions from the limma package (26). Gene coexpression networks were subsequently constructed to further explore the relationships among the prioritized risk genes, including TMEM190, CILP2, and FOXO3 (15). Briefly, the gene expression matrix of IVD was used to perform correlation analysis between each risk gene and all other genes. The genes were then ranked on the basis of their correlation indices. These correlation coefficients were then used for gene set enrichment analysis (GSEA) with pathway data from Reactome, which was performed via the clusterProfiler package and visualized via the Ridgeplot R package. This approach highlighted significant biological functions and pathways associated with each prioritized gene. Significant enrichment was determined on the basis of an adjusted p < 0.05, a normalized enrichment score (|NES|) > 1, and an FDR < 0.25.

PPI analysis

To investigate potential causal genes implicated in IDD, we employed the STRING database to perform extensive network analysis. The 88 proteins associated with IDD in both TWAS and PWAS were analyzed (Supplementary Data 1). We only reserved connections with an interaction score greater than 0.4.

Analysis of druggable genes and known drugs

To explore the druggability of potential causal genes of IDD, we conducted druggable gene and known drug analyses. A previous study categorized druggable genes into three tiers (27). We categorized the genes identified via TWAS and eQTL colocalization analyses, as well as genes identified via PWAS (ARIC and deCODE) and pQTL colocalization analyses, and further searched for updated information on the drugs targeting the identified putative causal proteins in the Open Target Platform,2 which is a comprehensive tool that promotes drug target identification via the integration of multiple databases.

Human IVD tissue collection

IVD tissue samples were collected from 6 patients undergoing spinal fusion surgery with the ethics approval of the Affiliated Hospital of Xuzhou Medical University (XYFY2023-KL337-01). The IDD cases were classified according to Pfirrmann’s method (28). Samples falling within the I–II classification were labeled as controls, while falling within the III–V classification were designated as severe IDD samples (29). Patient information is provided in Supplementary Table 1. Informed consent was obtained from all patients.

Western blot analysis

Total protein was extracted from IVD using a complete cell lysis buffer and quantified with the BCA protein assay kit (Beyotime, China). Protein samples were separated via sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) and transferred onto 0.2 μm PVDF membranes (Sigma-Aldrich, United States). The membranes were blocked with a 5% skim milk solution at room temperature and then incubated overnight at 4°C with primary antibodies specific to TMEM190 (1:500; Invitrogen, PA5-70986), CILP2 (1:500; Proteintech, 11813-1-AP), and FOXO3 (1:1,000; Biotime Biotechnology, AF609-1). Following three washes with Tris-buffered saline containing Tween 20 (TBST), the membranes were treated with secondary antibodies at room temperature. Immunoblotting was detected using the UVP ChemiDoc-It Imaging System (UVP, CA, United States) with an enhanced chemiluminescence detection kit (Thermo Fisher; 34,580) applied to the membranes. β-actin served as the loading control, and each blot was analyzed for integrated density using ImageJ software.

Animal experiments

All experiments were reviewed and approved by the committee of the Institutional Animal Care and Use Committee at Nanjing Drum Tower Hospital (approval number: DWSY-25005637). The methods were described as before (30). Briefly, after anesthesia, 31-G needle was poked into the 8-week-old male C57BL/6 IVD at 90° vertically, rotated 360° and held for 1 min. Sham operation was also performed. These operations were performed on coccygeal IVD. Immediately after the puncture, 2 ul shRNA targeted Cilp2 gene or negative control encapsulated with the GV112 vector (Shanghai Genechem Co., Ltd.) were injected into the IVD. Four weeks after acupuncture and shRNA injection, the caudal IVD of mice were examined by MRI. Degeneration grade of mice IVD was calculated as described before (28). IVD tissues were then collected for the test of Cilp2 protein levels and histological staining.

Histological staining and analysis

The harvested IVD tissue was immersed in a 4% paraformaldehyde solution (Solarbio, China) for 48 h to preserve its morphology. It was then decalcified in a 10% EDTA solution (pH 7.2–7.4) for 2 weeks, with daily changes of the solution. The tissue underwent dehydration through a series of graded ethanol baths and was subsequently cleared using an environmentally friendly clearing agent (Solarbio, China). After clearing, the tissue was embedded in paraffin wax and sectioned into 5 μm-thick slices using a microtome. Hematoxylin and Eosin (H & E) staining were applied according to the instructions (Solarbio, China) to observe the histological morphology of the IVD. Finally, histological scoring of the H & E samples was conducted following methodologies described before (31). Simply, the stained results were evaluated from two perspectives, including annulus fibrosus and nucleus pulposus. Each category was assigned a score ranging from 0 to 3, yielding a cumulative score between 0 and 6. Higher scoring levels indicated greater degrees of degeneration.

Statistical analyses

A chi-square test was conducted to assess the difference in the number of colocalizing genes between IDD-associated and non-associated genes identified by the TWAS. Furthermore, two-tailed t-test or one-way ANOVA test were employed to assess: differential expression of TMEM190, CILP2, and FOXO3 in degenerative compared to non-degenerative IVD using GEO dataset, protein levels of these genes in clinical and animal specimens, and degeneration grade, as well as H & E scores in a murine IDD model. Statistical significance (p < 0.05) is denoted by asterisks (*) in figures.

Results

Identification of genes associated with IDD by TWAS

As shown in the pipeline (Figure 1), we initiated our analysis by performing a cross-tissue TWAS based on data from reference gene expression predictions from GTEx and the largest intervertebral disc disorder GWAS conducted to date. TWAS analysis of 17,342 protein-coding genes identified 556 significantly IDD-associated genes (FDR-adjusted p < 0.05). The top 10 genes most strongly correlated with IDD were CHST3, SOX5, SPOCK2, SMAD3, FGFR3, C6orf106, GFPT1, TWIST1, ASCC1, IGFBP3 (Figure 2 and Supplementary Data 2). Among the 556 genes, 20 genes, including SMAD3, mapped to nearest genes at the previously reported GWAS susceptibility loci. Our analysis also suggested the novel associations of the remaining 536 genes with IDD risk.

FIGURE 2
Manhattan plot displaying genetic association data across chromosomes one to twenty-two. The x-axis represents chromosomes, and the y-axis shows the negative logarithm of the p-value. Red and blue labels indicate specific genes above and below statistical significance thresholds, respectively.

Figure 2. Manhattan plot illustrating TWAS gene associations. Each dot represents a gene plotted according to its genomic position (x-axis) and the significance of the association, measured as the −log10(FDR-adjusted p-value) (y-axis). Highlighted points with corresponding gene labels indicate genes meeting stringent colocalization criteria: FDR-adjusted p < 0.05 and colocalization max-GLCP > 0.5. The color of the highlighted points indicates the directionality of the genetic effect: red represents positive z-mean values (z-mean > 0) and blue for negative z-mean values (z-mean < 0).

We then conducted enrichment analysis of the significant findings using the KEGG, Reactome, and GO databases. Notably, the top enriched pathways of TWAS included extracellular matrix organization, cellular senescence, and skeletal system and connective tissue development, all of which established mechanisms in IDD pathogenesis. Other prominent enriched terms included calcineurin activates NFAT, glycosphingolipid biosynthesis, aspartate and asparagine metabolism, cartilage development, response to transforming growth factor-beta, chondrocyte differentiation, regulation of lipid kinase activity, and signal transduction pathways (Figures 3A, B and Supplementary Figure 1).

FIGURE 3
Charts displaying pathway enrichment analysis results, organized into six panels labeled A to F. Each panel shows pathways on the y-axis and the negative log of p-values on the x-axis, with red indicating high significance. Panels A and C highlight Reactome Enrichment for TWAS and ARIC, respectively; Panels B and D focus on KEGG Enrichment for TWAS and ARIC; Panels E and F show Reactome and KEGG Enrichment for deCODE. Pathways include extracellular matrix organization, sphingolipid metabolism, and bladder cancer, among others.

Figure 3. Reactome and KEGG enrichment analyses of the potential associated genes of IDD identified with TWAS and PWAS. Reactome and KEGG enrichment analysis of the potential associated genes of IDD identified with TWAS (A,B), with PWAS using ARIC dataset (C,D), and with PWAS using deCODE dataset (E,F). Each line represents a pathway with significance defined by an FDR-adjusted p < 0.05. The color intensity represents statistical significance. The dot size corresponds to the gene ratio, which is defined as the number of genes of a pathway to the total number of genes analyzed.

Colocalization between IDD risk loci and eQTLs

To determine whether TWAS-identified associations with IDD are driven by shared causal variants, we performed eQTL colocalization analysis using fastENLOC across 49 GTEx tissues for all protein-coding genes. This analysis identified 146 genes with strong colocalization evidence (max-GLCP > 0.5), among which 104 genes were TWAS- prioritized (Figure 2 and Supplementary Data 3, 4). TWAS-prioritized genes showed significant enrichment for colocalization signals compared to a matched background set (χ2 = 2195.6, fold-enrichment = 91.7, p < 0.001) (Supplementary Table 2), highlighting the specificity and robustness of our findings.

Identification of plasma proteins associated with IDD by PWAS

To identify proteins potentially associated with IDD for further validation, we conducted two independent PWAS by integrating GWAS summary statistics with human plasma proteomic data from the ARIC consortium and deCODE Genetics. The ARIC-based PWAS identified 494 significant associations, and the deCODE-based PWAS yielded 523 associations (Figures 4A, B; Supplementary Data 5, 6) (p < 0.05). Among these, 153 proteins were consistently associated with IDD across both datasets (Supplementary Data 7). Among them, six proteins, TMEM190, CILP2, FOXO3, SPON2, GALNT3, and NUF1, were additionally supported by TWAS and eQTL colocalization analyses, further reinforcing their relevance to IDD.

FIGURE 4
Manhattan plots displaying genetic associations. Panel A shows multiple significant peaks across various chromosomes with highlighted genes such as SERPINA1 and COMP. Panel B presents fewer significant peaks, with some genes like PLEKHA1 and RSPO2 marked. Data points represent -log10(p-values) plotted against chromosome numbers.

Figure 4. Manhattan plot illustrating PWAS protein associations on the basis of ARIC and deCODE data. Manhattan plot for the ARIC based (A) and deCODE based (B) PWAS of IDD. Each dot represents a protein plotted according to its genomic position (x-axis), and the significance of the association was measured as the −log10(p-value) (y-axis). Highlighted points and their protein labels indicate proteins meeting stringent colocalization criteria: p < 0.05 and colocalization PPH4 > 0.5. The color of the highlighted points indicates the directionality of the genetic effect: red for positive beta values (beta > 0) and blue for negative beta values (beta < 0).

For the PWAS results from the ARIC cohort, the most significant pathways were collagen formation, extracellular matrix organization, and glycosphingolipid/sphingolipid metabolism (Figures 3C, D and Supplementary Figure 2). Similarly, the deCODE PWAS analysis revealed glycosphingolipid/sphingolipid metabolism, regulation of actin cytoskeleton, extra-nuclear estrogen signaling, and signaling pathways associated with GPER1, NOTCH1, PI3K-Akt, and Hedgehog (Figures 3E, F and Supplementary Figure 3). The high degree of consistency between the TWAS and PWAS results across both datasets confirms the robustness of these findings.

Colocalization between IDD risk loci and pQTLs

To provide causal evidence for IDD-associated proteins, we performed pQTL colocalization analyses. In the ARIC dataset, 26 proteins demonstrated strong colocalization signals with IDD risk loci (PPH4 > 0.5; Supplementary Data 8). Among these, 22 proteins were also significantly associated with IDD in the PWAS (Figure 4A; Supplementary Data 9). In the deCODE dataset, pQTL colocalization identified 24 significant proteins (PPH4 > 0.5; Supplementary Data 10), of which 16 proteins showed consistent PWAS associations with IDD (Figure 4B; Supplementary Data 11). Totally, 10 proteins were identified as causal proteins via ARIC based and deCODE based PWAS and their respective colocalization analyses (Supplementary Data 12). Among these, TMEM190, CILP2, and FOXO3 were also supported by TWAS and eQTL colocalization analysis.

Collectively, TMEM190, CILP2, and FOXO3 emerged as proteins with strong causal evidence for IDD, supported across multiple omics layers including TWAS, two independent PWAS datasets, and both eQTL and pQTL colocalizations (Table 1).

TABLE 1
www.frontiersin.org

Table 1. Summary of three potential causal genes of IDD indicated by TWAS, PWAS and colocalization analyses.

Evaluation of the expression levels of genes/proteins identified by TWAS/PWAS

To identify differentially expressed genes (DEGs) in degenerative intervertebral discs, we analyzed mRNA expression profiles from human disc tissues using the microarray dataset GSE56081. Analysis of the GSE56081 dataset encompassing 13,170 genes captured 537 (96.6%), 395 (80.0%), and 404 (77.2%) of IDD-associated genes/proteins identified by TWAS, ARIC based, and deCODE based PWAS, respectively. Among these genes, 2,877 were significantly upregulated and 3,140 were downregulated in degenerated discs (Figure 5A; Supplementary Data 13). We observed 189 genes overlapped with TWAS-prioritized candidates (Supplementary Data 14) and 53 overlapped with proteins identified in both PWAS analyses (Supplementary Data 15).

FIGURE 5
Panel A shows a volcano plot highlighting genes TMEM190, CILP2, and FOXO3, with data points colored to indicate downregulated and upregulated groups. Panel B contains bar graphs comparing the expression levels of TMEM190, CILP2, and FOXO3 between control and IDD groups, with asterisk indicating significance. Panel C features normalized enrichment score (NES) plots for TMEM190, CILP2, and FOXO3, illustrating pathway involvement in biological processes according to Reactome Gene Set Enrichment Analysis (GSEA).

Figure 5. DEGs and enrichment analysis of the three potential causal genes of IDD. (A) Volcano plot of the DEG analysis. Each dot represents a gene plotted according to the significance of the association measured as the −log10 (FDR-adjusted p-value) (y-axis). The colors of the points are as follows: red for upregulated genes, blue for downregulated genes, and gray for non-differentially expressed genes (FDR-adjusted p ≥ 0.05). (B) Expression levels of potential causal genes for IDD in degenerative vs. non-degenerative groups, based on microassy data from clinical samples n = 5. Two-tailed t-test. *p < 0.05. (C) Significantly enriched pathways for the three potential causal genes (TMEM190, CILP2 and FOXO3) as determined by GSEA. Each line represents a pathway with significance defined by an FDR-adjusted p < 0.05. Yellow indicates upregulation (NES > 0), while green indicates downregulation (NES < 0).

In particular, all three genes (TMEM190, CILP2, and FOXO3) priorized by multiple-omics analyses were found to be significantly differentially expressed in degenerative intervertebral discs (p < 0.05) (Figures 5A, B and Table 1). Specifically, TMEM190 and CILP2 were upregulated in IDD samples, whereas FOXO3 was downregulated compared to the control group. Notably, CILP2 showed the most pronounced change, exhibiting a 1.5-fold increase in expression in degenerated discs relative to controls.

Functional annotation of TMEM190, CILP2, and FOXO3

The three potential causal genes were analyzed within the framework of GSEA to investigate their functions by coexpression analysis (Figure 5C). The GSEA results revealed that all of the three genes involved in the expression and translation of olfactory receptors, sensory perception and SRP-dependent cotranslational protein targeting to the membrane pathways.

To elucidate the interactions among the candidate genes associated with IDD (TMEM190, CILP2, and FOXO3), we performed PPI analysis involving 88 proteins associated with IDD identified in TWAS and PWAS (Supplementary Data 1). There were 27 genes whose connections had interaction scores greater than 0.4. Notably, TMEM190 did not interact with either CILP2 or FOXO3. Although no direct interactions were observed between CILP2 and FOXO3, several core proteins—SMAD3, COMP, IGFBP3, IGF1R, COL10A1, RUNX3, and PTK2—were identified as mediators of the interaction between CILP2 and FOXO3 (Figure 6).

FIGURE 6
Diagram of interconnected nodes representing gene interactions. Each node is labeled with gene names such as IGF1R, FOXO3, and SMAD3. Lines indicate interactions, highlighting key connections among the genes. The layout varies with some nodes having multiple connections.

Figure 6. PPI network of the significant proteins associated with IDD identified via TWAS and PWAS analyses. Each dot represents a protein. Lines denote physical PPIs.

Druggability of the identified genes and proteins

Among the genes identified through TWAS and eQTL analyses, 33 protein-coding genes were classified within the druggable genome: including 16 in tier 1, 8 in tier 2, and 9 in tier 3 (Supplementary Data 16). By searching the Open Target Platform, we identified several approved or investigational drugs targeting risk genes indicated by TWAS and eQTL analyses, including FGFR3, TGFA, CD79B, PDE3A, NQO1, AGER, ITGA2, CDK4, COL27A1, and PTK2 (Supplementary Data 17). Additionally, of the proteins identified through PWAS and pQTL analyses via ARIC or deCODE, 14 protein-coding genes were classified within the druggable genome: including 4 in tier 1 and 10 in tier 3 (Supplementary Data 18). We further identified several approved or investigational drugs targeting risk genes of IDD indicated by PWAS and pQTL analyses, namely PLG and PTHLH (Supplementary Data 19). Among the three potential causal genes, only CILP2 was druggable, while TMEM190 and FOXO3 were not in the druggable genome.

Validation studies for potential causal genes of IDD with clinical samples and animal model

To explore the roles of TMEM190, CILP2, and FOXO3 in IDD, we assessed their expression in human IVD specimens from mild degeneration (Grades I and II) and severe degeneration (Grades III, IV, and V). Western blot results showed increased expression of TMEM190 and CILP2 with concurrent decreased expression of FOXO3 in severely degenerated IVD tissues (Figures 7A, B), aligning with findings in DEGs analysis.

FIGURE 7
Side-by-side images with protein expression data. Panel A shows Western blots of CILP2, TMEM190, FOXO3, and β-ACTIN for control and severe IDD samples. Panel B presents bar graphs comparing relative protein expression levels for CILP2, TMEM190, and FOXO3, indicating significant increases in severe IDD compared to control.

Figure 7. Validation of risk genes in clinical IVD tissues. (A,B) Western-blot analysis of TMEM190, CILP2, and FOXO3 in control and severe IDD patients. n = 3. Data represent the mean (SD). Two tailed t-test. *p < 0.05.

Given the combined evidence supporting CILP2, and its classification as a druggable target, we further investigated its role in IDD using a needle-induced IVD degeneration mouse model. As shown in Figures 8A, B, the protein level of CILP2 in the IVD was up-regulated after puncture treatment. However, this level decreased significantly following shRNA transfection. MRI examinations revealed that the grade score of the IVD was significantly increased after puncture treatment, while it decreased markedly with CILP2 knockdown (KD) (Figures 8C,D). H & E staining showed the degenerated progression was alleviated following the down-regulation of CILP2 (Figures 8E, F). These results suggest that down-regulation of CILP2 reduces the susceptibility to intervertebral disc degeneration progression in the mouse model of IDD.

FIGURE 8
A series of panels depicting scientific data and analysis: A) Western blot showing expression of CILP2 and β-ACTIN across three conditions: Sham, Puncture, and KD. B) Bar graph illustrating relative protein expression, with significant differences marked with asterisks and “ns” indicating no significance between Sham and Puncture. C) MRI images comparing control and KD conditions for puncture and sham treatments. D) Bar graph of Pfirrmann grade across conditions, with significant differences indicated. E) Histological images of tissue under Sham, Puncture, and Puncture + KD conditions. F) Bar graph showing histological scores with significant differences marked.

Figure 8. Down-regulation of CILP2 alleviated IDD progression in mouse model of IDD. (A,B) Western-blot analysis of CILP2 in Sham, puncture (Punc) and knocking down (KD) IVD tissues in mouse model of IDD. n = 3. (C,D) Magnetic resonance imaging (MRI) and Pfirrmann grades of IVD in mice treated as in (A). n = 3. (E,F) Hematoxylin and Eosin (H&E) staining and histological score assessment of IVD in mice treated as in (A). n = 3. Scale bar = 200 μm. Data represent the mean (SD). One-way ANOVA. *p < 0.05.

Discussion

To the best of our knowledge, this study is the first to employ multidimensional multi-omics data, including high-throughput genomics, whole-body transcriptomics, plasma proteomics and intervertebral disc transcriptomics, to investigate potential risk genes for IDD. Our integrative approach presented 104 TWAS-identified genes and 10 PWAS-identified proteins with IDD based on converging evidence supported by eQTL/pQTL colocalization analyses. These genes/proteins were enriched for key regulators of disc pathology, such as glycosphingolipid/sphingolipid metabolism. Three potential causal genes, TMEM190, CILP2, and FOXO3, were consistently supported by TWAS, two independent PWAS and colocalization analyses. These three genes were dysregulated in degenerated human discs, with CILP2 further classified as druggable. We also validated the role of these causal genes, TMEM190, CILP2 and FOXO3 with clinical samples, as well as the role of CILP2 with animal model in IDD.

Glycosphingolipid/sphingolipid metabolism consistently emerged as a key pathway across all enrichment analyses of the identified associations. Both the TWAS and the PWAS results from the ARIC and deCODE cohorts strongly highlight this pathway as a critical factor in the pathogenesis of IDD. Sphingolipids, including ceramide and sphingosine-1-phosphate, constitute a major class of lipids found in all eukaryotic cells. Sphingolipids regulate a wide range of biological processes, including inflammation, mitochondrial function and apoptosis (3234). The metabolic processes involved in sphingolipid biosynthesis and regulation were significantly enriched, underscoring their potential role in IDD. This pathway’s involvement in inflammation and apoptosis suggests that targeting sphingolipids synthesis could serve as a promising therapeutic strategy to alleviate disc disorders and associated pain.

Among the identified genes (TMEM190, CILP2, and FOXO3), FOXO3 has been previously investigated in the context of IDD. FOXO3, a member of the forkhead box O transcription factor family, is known to regulate critical cellular processes, including the cell cycle, apoptosis, and metabolism, and is implicated in age-related diseases (35, 36). FOXO3 has been linked to IDD in numerous studies, where it functions as a mediator regulating the role of specific genes in the disease, such as YTHDF2 and P300 (37, 38). Furthermore, FOXO3 is involved in the molecular mechanisms of potential therapeutic agents for IDD, such as stem cell-derived exosomes and procyanidin C1, primarily by regulating oxidative stress (39, 40), which reinforces the potential of FOXO3 as therapeutic interventions for this condition. The current study found that FOXO3 was down-regulated in severely degenerated disc tissues, which aligns with previous findings (38, 41). Our study provides additional evidence for dysregulated FOXO3 expression in IDD. However, the contribution of FOXO3 to IDD still requires further exploration.

The roles of the other two identified genes, TMEM190 and CILP2, in IDD are less well characterized. TMEM190 is located on chromosome 7 and contains five exons, which encode a small single-pass transmembrane protein (42). Small single-pass transmembrane proteins may be associated with mitochondrial oxidative phosphorylation (43), which has been linked to IDD. Additionally, TMEM190 may contribute to chondrocyte dedifferentiation (44). Given that cartilage endplate remodeling and altered chondrocyte subsets play key roles in IDD progression (45, 46), TMEM190 may involve in IDD pathogenesis. CILP2, a member of the cartilage intermediate layer protein family, encodes a matricellular protein predominantly expressed in cartilage cells but also in various other tissues (47). Quantitative proteomic analysis and immunohistochemistry have demonstrated increased CILP2 levels in degenerated human intervertebral discs (48, 49). In current study, we found CILP2 was up-regulated in severe IDD tissues, which provides additional evidence that CILP2 play a role in IDD. Importantly, the inhibition of Cilp2 has been shown to improve mitochondrial dysfunction in sarcopenia via the WNT signaling pathway (47). Given the established roles of mitochondrial dysfunction and WNT signaling in IDD (50, 51), CILP2 is likely to play an important role in this condition. In the current study, we found that down-regulation of CILP2 alleviated IDD progression in mouse model of IDD. Our results provide more direct evidence for the role of CILP2 in the progression of IDD. Our multi-omics investigation and validation study with clinical samples and animal model offer evidence supporting the role of CILP2 as a disease-causing gene and therapeutic target in IDD. However, the functional mechanism of CILP2 in IDD was not explored in current study. It has been reported that CILP2 affect sarcopenia and hypertrophic scar by antagonizing Wnt signaling pathway (47), and reducing the ubiquitination of ACLY (52), respectively. Further research is warranted to elucidate the precise role of CILP2 in modulating IDD progression through these candidate signaling pathways.

While no direct interactions among the three genes (TMEM190, CILP2 and FOXO3) were reported, GSEA revealed their collective involvement in the olfactory signaling pathway, and sensory perception. Notably, olfactory stem cells have been shown to exhibit a chondrogenic phenotype, promoting IVD regeneration in a rat model of disc injury (53), which indicates these three genes collectively contributed to the pathology of IDD. In addition, pain is a significant symptom of IDD. It has been reported that anti-sensory nerve transmission significantly suppresses inflammatory pain markers (54). The involvement of these three genes in sensory perception indicates that they share a similar pathway for the contribution of pain to IDD. To further explore the correlations, we performed a PPI analysis using data from TWAS and PWAS. Although no direct interactions were found among the three genes, several mediators of interaction between CILP2 and FOXO3 were identified, including COMP and IGFBP3, which have been linked to IDD progression (55, 56). These findings suggest that CILP2 and FOXO3 may collaboratively influence IDD through these mediators.

Our study has several strengths. Primarily, it integrates both genomic and proteomic data to provide comprehensive insights into the complex biological systems underlying IDD. Additionally, the validation of potential causal genes through two independent PWAS analyses strengthens the reliability of our findings. Furthermore, the datasets utilized, comprising extensive human transcriptomes, proteomes, and IDD GWAS data, are among the largest and most complete to date, enhancing the robustness of the results. Finally, we validated the risk genes derived from public datasets using clinical samples and animal model, which enhances the reliability of our findings.

Some limitations must be acknowledged. First, the PWAS for IDD utilized human blood proteome data; however, plasma proteins serve as systemic biomarkers and may not fully capture disc-specific changes, potentially introducing bias. Future studies should examine the proteomes specific to human intervertebral discs. Secondly, the identification of susceptible genes from a European database, coupled with clinical validation using specimens from the Chinese population, introduces population heterogeneity that may limit the generalizability of the findings. Future cross-ethnic validation studies should be conducted to assess the robustness of these findings across diverse populations and ensure their applicability in broader clinical contexts. Additionally, the mechanisms by which identified risk genes and the relevance of their enriched pathways contribute to IDD remain unclear, and additional studies are needed to further evaluate their potential as therapeutic targets. Besides, only the mouse tail disc puncture model was employed, which is an acute injury model and may not adequately replicate the chronic, progressive nature of human IDD. Also, the current transcriptomic samples predominantly represent European populations, whereas the proteomic samples are from American populations, and expanding the diversity of these datasets may yield more accurate estimations and broader applicability. Finally, gene-environment interactions and assortative mating could influence genetic effects and contribute to variance in the analysis. Unfortunately, due to the limitations inherent in the current dataset and the scope of the study design, it is not feasible to adjust for these factors in this particular analysis. Nevertheless, the strength of our study lies in its innovative integration of multi-omics data, which positions it as one of the first efforts to identify and validate novel genetic risk factors for IDD in such a comprehensive manner.

Conclusion

In summary, we present an expanded resource of putatively causal genes associated with IDD, and highlight three novel potential causal genes (TMEM190, CILP2, and FOXO3). These findings provided a broad hint for further research on the potential mechanisms underlying IDD pathogenesis and highlight novel therapeutic targets for future investigations.

Data availability statement

All the data used in this study are publicly available without the need for special access. The specific sources of the data are as follows: TWAS Process: https://github.com/xqwen/fastenloc. IDD GWAS: https://pmc.ncbi.nlm.nih.gov/articles/PMC8810832. TWAS Prediction Model: https://predictdb.org/post/2021/07/21/gtex-v8-models-on-eqtl-and-sqtl/. fastENLOC Colocalization: Pre-computed GTEx multi-tissue eQTL annotations with hg38 position ID from https://github.com/xqwen/fastenloc. BLISS Software: Data calculations were performed using BLISS, which can be accessed at https://github.com/gcb-hub/BLISS. GEO database (GSE56081).

Ethics statement

The studies involving humans were approved by the Committee of the Affiliated Hospital of Xuzhou Medical University. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin. The animal study was approved by the committee of the Institutional Animal Care and Use Committee at Nanjing Drum Tower Hospital. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

ZZha: Writing – original draft. YC: Writing – review & editing. QW: Writing – review & editing. ZL: Writing – review & editing. BD: Writing – review & editing. CD: Writing – review & editing. ZZhu: Writing – review & editing.

Funding

The author(s) declare financial support was received for the research and/or publication of this article. This work was supported by the Talent Introduction Project of the Affiliated Hospital of Xuzhou Medical University (2024203033), the Science and Technology Project of the Affiliated Hospital of Xuzhou Medical University (2024ZL32), and the National Nature Science Foundation of China (82202700).

Acknowledgments

We thank all the participants and investigators of GWAS, transcriptomics and proteomics studies included in this study.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The authors declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmed.2025.1698050/full#supplementary-material

Supplementary Figure 1 | GO enrichment analysis of the potential associated genes of IDD identified with TWAS. Each line represents a pathway with significance defined by an FDR-adjusted p < 0.05. The color intensity represents statistical significance. The dot size corresponds to the gene ratio, which is defined as the number of genes of a pathway to the total number of genes analyzed.

Supplementary Figure 2 | GO enrichment analysis of the potential associated genes of IDD identified with PWAS via the ARIC dataset. Each line represents a pathway with significance defined by an FDR-adjusted p < 0.05. The color intensity represents statistical significance. The dot size corresponds to the gene ratio, which is defined as the number of genes of a pathway to the total number of genes analyzed.

Supplementary Figure 3 | GO enrichment analysis of the potential associated genes of IDD identified with PWAS via the deCODE dataset. Each line represents a pathway with significance defined by an FDR-adjusted p < 0.05. The color intensity represents statistical significance. The dot size corresponds to the gene ratio, which is defined as the number of genes of a pathway to the total number of genes analyzed.

Footnotes

References

1. Livshits G, Popham M, Malkin I, Sambrook P, Macgregor A, Spector T, et al. Lumbar disc degeneration and genetic factors are the main risk factors for low back pain in women: the UK Twin Spine study. Ann Rheum Dis. (2011) 70:1740–5. doi: 10.1136/ard.2010.137836

PubMed Abstract | Crossref Full Text | Google Scholar

2. Xu J, Shao T, Lou J, Zhang J, Xia C. Aging, cell senescence, the pathogenesis and targeted therapies of intervertebral disc degeneration. Front Pharmacol. (2023) 14:1172920. doi: 10.3389/fphar.2023.1172920

PubMed Abstract | Crossref Full Text | Google Scholar

3. Ou-Yang D, Kleck C, Ackert-Bicknell C. Genetics of intervertebral disc degeneration. Curr Osteoporos Rep. (2023) 21:56–64. doi: 10.1007/s11914-022-00769-0

PubMed Abstract | Crossref Full Text | Google Scholar

4. Vergroesen P, Kingma I, Emanuel K, Hoogendoorn R, Welting T, van Royen B, et al. Mechanics and biology in intervertebral disc degeneration: a vicious circle. Osteoarthritis Cartilage. (2015) 23:1057–70. doi: 10.1016/j.joca.2015.03.028

PubMed Abstract | Crossref Full Text | Google Scholar

5. Guo W, Li B, Zhao J, Li X, Wang L. Causal associations between modifiable risk factors and intervertebral disc degeneration. Spine J. (2024) 24:195–209. doi: 10.1016/j.spinee.2023.10.021

PubMed Abstract | Crossref Full Text | Google Scholar

6. Foster N, Anema J, Cherkin D, Chou R, Cohen S, Gross D, et al. Prevention and treatment of low back pain: evidence, challenges, and promising directions. Lancet. (2018) 391:2368–83. doi: 10.1016/S0140-6736(18)30489-6

PubMed Abstract | Crossref Full Text | Google Scholar

7. Bjornsdottir G, Stefansdottir L, Thorleifsson G, Sulem P, Norland K, Ferkingstad E, et al. Rare SLC13A1 variants associate with intervertebral disc disorder highlighting role of sulfate in disc pathology. Nat Commun. (2022) 13:634. doi: 10.1038/s41467-022-28167-1

PubMed Abstract | Crossref Full Text | Google Scholar

8. Suri P, Naeini M, Heagerty P, Freidin M, Smith I, Elgaeva E, et al. The association of lumbar intervertebral disc degeneration with low back pain is modified by underlying genetic propensity to pain. Spine J. (2025) 25:8–17. doi: 10.1016/j.spinee.2024.05.018

PubMed Abstract | Crossref Full Text | Google Scholar

9. Mai J, Lu M, Gao Q, Zeng J, Xiao J. Transcriptome-wide association studies: recent advances in methods, applications and available databases. Commun Biol. (2023) 6:899. doi: 10.1038/s42003-023-05279-y

PubMed Abstract | Crossref Full Text | Google Scholar

10. Nica A, Dermitzakis E. Expression quantitative trait loci: present and future. Philos Trans R Soc Lond B Biol Sci. (2013) 368:20120362. doi: 10.1098/rstb.2012.0362

PubMed Abstract | Crossref Full Text | Google Scholar

11. GTEx Consortium. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science. (2020) 369:1318–30. doi: 10.1126/science.aaz1776

PubMed Abstract | Crossref Full Text | Google Scholar

12. GTEx Consortium. Genetic effects on gene expression across human tissues. Nature. (2017) 550:204–13. doi: 10.1038/nature24277

PubMed Abstract | Crossref Full Text | Google Scholar

13. Zhang J, Dutta D, Köttgen A, Tin A, Schlosser P, Grams M, et al. Plasma proteome analyses in individuals of European and African ancestry identify cis-pQTLs and models for proteome-wide association studies. Nat Genet. (2022) 54:593–602. doi: 10.1038/s41588-022-01051-w

PubMed Abstract | Crossref Full Text | Google Scholar

14. Wingo A, Liu Y, Gerasimov E, Gockley J, Logsdon B, Duong D, et al. Integrating human brain proteomes with genome-wide association data implicates new proteins in Alzheimer’s disease pathogenesis. Nat Genet. (2021) 53:143–6. doi: 10.1038/s41588-020-00773-z

PubMed Abstract | Crossref Full Text | Google Scholar

15. Wu C, Liu H, Zuo Q, Jiang A, Wang C, Lv N, et al. Identifying novel risk genes in intracranial aneurysm by integrating human proteomes and genetics. Brain. (2024) 147:2817–25. doi: 10.1093/brain/awae111

PubMed Abstract | Crossref Full Text | Google Scholar

16. Dai Z, Wu Y, Huang H, Zheng H. Integrating brain proteomes and genetics to identify novel risk genes in chronic widespread musculoskeletal pain. Sci Rep. (2025) 15:21999. doi: 10.1038/s41598-025-04379-5

PubMed Abstract | Crossref Full Text | Google Scholar

17. Ferkingstad E, Sulem P, Atlason B, Sveinbjornsson G, Magnusson M, Styrmisdottir E, et al. Large-scale integration of the plasma proteome with genetics and disease. Nat Genet. (2021) 53:1712–21. doi: 10.1038/s41588-021-00978-w

PubMed Abstract | Crossref Full Text | Google Scholar

18. Barbeira A, Pividori M, Zheng J, Wheeler H, Nicolae D, Im H. Integrating predicted transcriptome from multiple tissues improves association detection. PLoS Genet. (2019) 15:e1007889. doi: 10.1371/journal.pgen.1007889

PubMed Abstract | Crossref Full Text | Google Scholar

19. Al-Barghouthi B, Rosenow W, Du K, Heo J, Maynard R, Mesner L, et al. Transcriptome-wide association study and eQTL colocalization identify potentially causal genes responsible for human bone mineral density GWAS associations. Elife. (2022) 11:e77285. doi: 10.7554/eLife.77285

PubMed Abstract | Crossref Full Text | Google Scholar

20. Yuan L, Su Y, Zhao J, Cho M, Wang D, Yuan L, et al. Investigating the shared genetic architecture between obesity and depression: a large-scale genomewide cross-trait analysis. Front Endocrinol. (2025) 16:1578944. doi: 10.3389/fendo.2025.1578944

PubMed Abstract | Crossref Full Text | Google Scholar

21. Giambartolomei C, Vukcevic D, Schadt E, Franke L, Hingorani A, Wallace C, et al. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genet. (2014) 10:e1004383. doi: 10.1371/journal.pgen.1004383

PubMed Abstract | Crossref Full Text | Google Scholar

22. Foley C, Staley J, Breen P, Sun B, Kirk P, Burgess S, et al. A fast and efficient colocalization algorithm for identifying shared genetic risk factors across multiple traits. Nat Commun. (2021) 12:764. doi: 10.1038/s41467-020-20885-8

PubMed Abstract | Crossref Full Text | Google Scholar

23. Zuber V, Grinberg N, Gill D, Manipur I, Slob E, Patel A, et al. Combining evidence from Mendelian randomization and colocalization: review and comparison of approaches. Am J Hum Genet. (2022) 109:767–82. doi: 10.1016/j.ajhg.2022.04.001

PubMed Abstract | Crossref Full Text | Google Scholar

24. Yu G, Wang L, Han Y, He Q. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | Crossref Full Text | Google Scholar

25. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. (2000) 28:27–30. doi: 10.1093/nar/28.1.27

PubMed Abstract | Crossref Full Text | Google Scholar

26. Ritchie M, Phipson B, Wu D, Hu Y, Law C, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

27. Finan C, Gaulton A, Kruger F, Lumbers R, Shah T, Engmann J, et al. The druggable genome and support for target identification and validation in drug development. Sci Transl Med. (2017) 9:eaag1166. doi: 10.1126/scitranslmed.aag1166

PubMed Abstract | Crossref Full Text | Google Scholar

28. Pfirrmann C, Metzdorf A, Zanetti M, Hodler J, Boos N. Magnetic resonance classification of lumbar intervertebral disc degeneration. Spine. (2001) 26:1873–8. doi: 10.1097/00007632-200109010-00011

PubMed Abstract | Crossref Full Text | Google Scholar

29. Chen S, Lei L, Li Z, Chen F, Huang Y, Jiang G, et al. Grem1 accelerates nucleus pulposus cell apoptosis and intervertebral disc degeneration by inhibiting TGF-β-mediated Smad2/3 phosphorylation. Exp Mol Med. (2022) 54:518–30. doi: 10.1038/s12276-022-00753-9

PubMed Abstract | Crossref Full Text | Google Scholar

30. Jiang H, Qin H, Yang Q, Huang L, Liang X, Wang C, et al. Effective delivery of miR-150-5p with nucleus pulposus cell-specific nanoparticles attenuates intervertebral disc degeneration. J Nanobiotechnology. (2024) 22:292. doi: 10.1186/s12951-024-02561-x

PubMed Abstract | Crossref Full Text | Google Scholar

31. Keorochana G, Johnson J, Taghavi C, Liao J, Lee K, Yoo J, et al. The effect of needle size inducing degeneration in the rat caudal disc: evaluation using radiograph, magnetic resonance imaging, histology, and immunohistochemistry. Spine J. (2010) 10:1014–23. doi: 10.1016/j.spinee.2010.08.013

PubMed Abstract | Crossref Full Text | Google Scholar

32. York A, Skadow M, Oh J, Qu R, Zhou Q, Hsieh W, et al. IL-10 constrains sphingolipid metabolism to limit inflammation. Nature. (2024) 627:628–35. doi: 10.1038/s41586-024-07098-5

PubMed Abstract | Crossref Full Text | Google Scholar

33. Hernández-Corbacho M, Salama M, Canals D, Senkal C, Obeid L. Sphingolipids in mitochondria. Biochim Biophys Acta Mol Cell Biol Lipids. (2017) 1862:56–68. doi: 10.1016/j.bbalip.2016.09.019

PubMed Abstract | Crossref Full Text | Google Scholar

34. Tirodkar T, Voelkel-Johnson C. Sphingolipids in apoptosis. Exp Oncol. (2012) 34:231–42.

Google Scholar

35. Eijkelenboom A, Burgering B. FOXOs: signalling integrators for homeostasis maintenance. Nat Rev Mol Cell Biol. (2013) 14:83–97. doi: 10.1038/nrm3507

PubMed Abstract | Crossref Full Text | Google Scholar

36. McIntyre RL, Liu Y, Hu M, Morris B, Willcox B, Donlon T, et al. Pharmaceutical and nutraceutical activation of FOXO3 for healthy longevity. Ageing Res Rev. (2022) 78:101621. doi: 10.1016/j.arr.2022.101621

PubMed Abstract | Crossref Full Text | Google Scholar

37. Wang F, Wang Y, Zhang S, Pu M, Zhou P. YTHDF2-dependent m6A modification of FOXO3 mRNA mediates TIMP1 expression and contributes to intervertebral disc degeneration following ROS stimulation. Cell Mol Life Sci. (2024) 81:477. doi: 10.1007/s00018-024-05503-w

PubMed Abstract | Crossref Full Text | Google Scholar

38. Hao Y, Ren Z, Yu L, Zhu G, Zhang P, Zhu J, et al. p300 arrests intervertebral disc degeneration by regulating the FOXO3/Sirt1/Wnt/β-catenin axis. Aging Cell. (2022) 21:e13677. doi: 10.1111/acel.13677

PubMed Abstract | Crossref Full Text | Google Scholar

39. Bian Z, Zhai Y, Zhang Y, Wang T, Li H, Ouyang J, et al. Senescent cartilage endplate stem cells-derived exosomes induce oxidative stress injury in nucleus pulposus cells and aggravate intervertebral disc degeneration by regulating FOXO3. Free Radic Biol Med. (2025) 233:39–54. doi: 10.1016/j.freeradbiomed.2025.03.027

PubMed Abstract | Crossref Full Text | Google Scholar

40. Hua W, Xie L, Dong C, Yang G, Chi S, Xu Z, et al. Procyanidin C1 ameliorates acidic pH stress induced nucleus pulposus degeneration through SIRT3/FOXO3-mediated mitochondrial dynamics. J Transl Med. (2024) 22:1071. doi: 10.1186/s12967-024-05805-4

PubMed Abstract | Crossref Full Text | Google Scholar

41. Xia P, Gao X, Li F, Shao L, Sun Y. Down-regulation of microRNA-30d alleviates intervertebral disc degeneration through the promotion of FOXO3 and suppression of CXCL10. Calcif Tissue Int. (2021) 108:252–64. doi: 10.1007/s00223-020-00760-w

PubMed Abstract | Crossref Full Text | Google Scholar

42. Nishimura H, Gupta S, Myles D, Primakoff P. Characterization of mouse sperm TMEM190, a small transmembrane protein with the trefoil domain: evidence for co-localization with IZUMO1 and complex formation with other sperm proteins. Reproduction. (2011) 141:437–51. doi: 10.1530/REP-10-0391

PubMed Abstract | Crossref Full Text | Google Scholar

43. Zickermann V, Angerer H, Ding M, Nübel E, Brandt U. Small single transmembrane domain (STMD) proteins organize the hydrophobic subunits of large membrane protein complexes. FEBS Lett. (2010) 584:2516–25. doi: 10.1016/j.febslet.2010.04.021

PubMed Abstract | Crossref Full Text | Google Scholar

44. Lindberg ED, Kaya S, Jamali A, Alliston T, O’Connell G. Effect of passaging on bovine chondrocyte gene expression and engineered cartilage production. Tissue Eng Part A. (2024) 30:512–24. doi: 10.1089/ten.TEA.2023.0349

PubMed Abstract | Crossref Full Text | Google Scholar

45. Li H, Tang Y, Liu Z, Chen K, Zhang K, Hu S, et al. Lumbar instability remodels cartilage endplate to induce intervertebral disc degeneration by recruiting osteoclasts via Hippo-CCL3 signaling. Bone Res. (2024) 12:34. doi: 10.1038/s41413-024-00331-x

PubMed Abstract | Crossref Full Text | Google Scholar

46. Zhang Y, Han S, Kong M, Tu Q, Zhang L, Ma X. Single-cell RNA-seq analysis identifies unique chondrocyte subsets and reveals involvement of ferroptosis in human intervertebral disc degeneration. Osteoarthritis Cartilage. (2021) 29:1324–34. doi: 10.1016/j.joca.2021.06.010

PubMed Abstract | Crossref Full Text | Google Scholar

47. Deng Z, Song C, Chen L, Zhang R, Yang L, Zhang P, et al. Inhibition of CILP2 improves glucose metabolism and mitochondrial dysfunction in sarcopenia via the Wnt signalling pathway. J Cachexia Sarcopenia Muscle. (2024) 15:2544–58. doi: 10.1002/jcsm.13597

PubMed Abstract | Crossref Full Text | Google Scholar

48. Yee A, Lam M, Tam V, Chan W, Chu I, Cheah K, et al. Fibrotic-like changes in degenerate human intervertebral discs revealed by quantitative proteomic analysis. Osteoarthritis Cartilage. (2016) 24:503–13. doi: 10.1016/j.joca.2015.09.020

PubMed Abstract | Crossref Full Text | Google Scholar

49. Koiv K, Aunapuu M, Torga T, Rätsep T, Bakhoff K, Arend A, et al. CILP-2 expression in the intervertebral discs of patients with lumbar radiculopathy. BMC Musculoskelet Disord. (2024) 25:882. doi: 10.1186/s12891-024-07996-9

PubMed Abstract | Crossref Full Text | Google Scholar

50. Song Y, Liang H, Li G, Ma L, Zhu D, Zhang W, et al. The NLRX1-SLC39A7 complex orchestrates mitochondrial dynamics and mitophagy to rejuvenate intervertebral disc by modulating mitochondrial Zn(2+) trafficking. Autophagy. (2024) 20:809–29. doi: 10.1080/15548627.2023.2274205

PubMed Abstract | Crossref Full Text | Google Scholar

51. Dong W, Liu J, Lv Y, Wang F, Liu T, Sun S, et al. miR-640 aggravates intervertebral disc degeneration via NF-kappaB and WNT signalling pathway. Cell Prolif. (2019) 52:e12664. doi: 10.1111/cpr.12664

PubMed Abstract | Crossref Full Text | Google Scholar

52. Wang J, Du J, Wang Y, Song Y, Wu J, Wang T, et al. CILP2 promotes hypertrophic scar through Snail acetylation by interaction with ACLY. Biochim Biophys Acta Mol Basis Dis. (2024) 1870:167202. doi: 10.1016/j.bbadis.2024.167202

PubMed Abstract | Crossref Full Text | Google Scholar

53. Murrell W, Sanford E, Anderberg L, Cavanagh B, Mackay-Sim A. Olfactory stem cells can be induced to express chondrogenic phenotype in a rat intervertebral disc injury model. Spine J. (2009) 9:585–94. doi: 10.1016/j.spinee.2009.02.011

PubMed Abstract | Crossref Full Text | Google Scholar

54. Nojima D, Inage K, Sakuma Y, Sato J, Orita S, Yamauchi K, et al. Efficacy of anti-NaV1.7 antibody on the sensory nervous system in a rat model of lumbar intervertebral disc injury. Yonsei Med J. (2016) 57:748–53. doi: 10.3349/ymj.2016.57.3.748

PubMed Abstract | Crossref Full Text | Google Scholar

55. Ding JY, Yan X, Zhang R, Zhang H, Kang L, Jia C, et al. Diagnostic value of serum COMP and ADAMTS7 for intervertebral disc degeneration. Eur J Med Res. (2024) 29:196. doi: 10.1186/s40001-024-01784-w

PubMed Abstract | Crossref Full Text | Google Scholar

56. Kazezian Z, Li Z, Alini M, Grad S, Pandit A. Injectable hyaluronic acid down-regulates interferon signaling molecules, IGFBP3 and IFIT3 in the bovine intervertebral disc. Acta Biomater. (2017) 52:118–29. doi: 10.1016/j.actbio.2016.12.029

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: intervertebral disc disorder, transcriptome-wide association study, proteome-wide association study, genome-wide association study, validation study

Citation: Zhang Z, Chen Y, Wang Q, Li Z, Dai B, Dong C and Zhu Z (2025) Identification and validation of novel risk genes for intervertebral disc disorder by integrating large-scale multi-omics analyses and experimental studies. Front. Med. 12:1698050. doi: 10.3389/fmed.2025.1698050

Received: 03 September 2025; Accepted: 27 October 2025;
Published: 12 November 2025.

Edited by:

HaiHui Huang, Shaoguan University, China

Reviewed by:

Xuan Song, Harvard University, United States
Tuo Han, The Second Affiliated Hospital of Xi’an Jiaotong University, China

Copyright © 2025 Zhang, Chen, Wang, Li, Dai, Dong and Zhu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Zhengya Zhu, eHlmeXp6eTI4MUAxMjYuY29t

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.