Multi-Omics Bioinformatic Analyses Linking LOXL4 with Spondylolisthesis

Introduction

Spondylolisthesis refers to the anterior or posterior slippage of one vertebra relative to an adjacent vertebra, typically occurring in the lumbar region.¹ The etiology of this condition is complex, and involves factors such as intervertebral disc degeneration, facet joint instability, laxity of surrounding ligaments, congenital developmental anomalies, and trauma. The clinical manifestations primarily include lower back pain and intermittent claudication due to nerve compression; in severe cases, it can lead to paralysis.^2,3 The treatment strategy is contingent upon the degree of slippage, nerve compression, and the patient’s overall health status. Conservative treatments, which include physical therapy, medication, and lifestyle modifications, are aimed at alleviating symptoms and enhancing quality of life; however, they do not fundamentally address the mechanical instability and neurological dysfunction associated with spondylolisthesis.⁴ Procedures such as spinal canal decompression and spinal fusion are indicated for patients with severe symptoms who do not respond to conservative treatment. An article published by Seip A et al⁵ in JAMA highlights the significant controversy surrounding the decision to perform fusion based solely on decompression, moreover surgical interventions also entail certain risks and limitations, particularly concerning postoperative complications and functional recovery.

Previous studies have identified various risk factors for spondylolisthesis, including increasing age, high body mass index (BMI), intervertebral disc degeneration, facet joint degeneration, and reduced paraspinal muscle content.^6,7 Among these, disc degeneration is widely regarded as a key prerequisite for spondylolisthesis.⁸ At the cellular level, disc degeneration is thought to be initiated by an imbalance between extracellular matrix synthesis and degradation, along with a decline in cellular function.^9,10 Local inflammatory infiltration and oxidative stress accelerate matrix degradation, primarily through the activation of matrix metalloproteinases (MMPs) triggered by pro-inflammatory factors such as TNF-α and IL-1β. This cascade induces apoptosis of nucleus pulposus, leads to proteoglycan loss, and promotes collagen accumulation within the nucleus pulposus, ultimately destabilizing the spine and contributing to the development of spondylolisthesis.¹¹ In turn, the progressive disc degeneration often accompanied by facet joint degeneration compromises spinal stability and can precipitate spondylolisthesis. Conversely, once a vertebral slip occurs, the resulting aberrant load distribution imposes high shear stress on adjacent disc, which accelerates further disc degeneration in a self-perpetuating vicious cycle of instability and degeneration.⁷

In addition to these biomechanical factors, genetic predisposition plays a crucial role in the pathogenesis of spondylolisthesis. A systematic review by Mayer et al,¹² highlighted several candidate gene polymorphisms associated with intervertebral disc degeneration and related spinal disorders, such as ACAN (aggrecan), collagen type IX (COL9), GDF5 (growth differentiation factor 5), and ASPN (asporin), which collectively influence the pathogenesis of spinal degenerative diseases. Furthermore, genetic studies have also demonstrated that mutation in genes including CDMP-1, GFPT1, NFU1, AAK1, and LOC genes heighten the risk of spondylolisthesis.^13–17 As such, spondylolisthesis is not merely a biomechanical disorder but is closely linked to genetic and molecular factors. Genetic research plays a pivotal role in uncovering the disease’s underlying mechanisms, thus advancing precision medicine.

Summary-data-based Mendelian Randomization (SMR) analysis is a novel multi-omics data integration method proposed by Zhu et al,¹⁸ based on the principle of two-sample Mendelian randomization. This approach employs genetic variants as instrumental variables, making use of genome-wide association study (GWAS) summary data to investigate the causal links between exposure factors and disease outcomes. SMR analysis is grounded in three primary assumptions: a significant correlation exists between the instrumental variables Single Nucleotide Polymorphisms (SNPs) and the exposure factors, which ensures that genetic variations adequately represent these exposure factors; the instrumental variables SNPs do not have a direct connection to potential confounding factors, thus preventing bias in causal inferences; the effect of the instrumental variables SNPs on the outcome is mediated exclusively through the exposure factors and does not impact the outcome directly.¹⁹ In contrast to single MR analysis, multi-omics MR analysis amalgamates data from various omics levels including genomics, transcriptomics, and epigenetics, allowing for a more thorough examination of the diverse effects of genetic factors on diseases and thereby improving the accuracy and explanatory capacity of the analysis. SMR analysis depends on GWAS summary statistics, providing high privacy protection levels and benefits in sharing data. Consequently, it possesses significant application value in public databases and aids in the exploration of the genetic mechanisms that underpin complex diseases.²⁰

This study utilizes a multi-omics strategy to clarify the molecular foundations and clinical significance of spondylolisthesis development. By using genomic datasets related to spondylolisthesis from public databases, we applied stringent bioinformatics techniques, such as SMR analysis and Bayesian colocalization analysis, to pinpoint high-confidence causal genes linked to the pathology of vertebral slippage. It is important to note that the eQTL and mQTL datasets are derived from peripheral blood samples, not from intervertebral disc or other spinal tissues in our SMR analysis. This choice was necessitated by the lack of any publicly available eQTL or mQTL resources specific to human spinal or disc tissue.^21,22 Indeed, current QTL studies are often biased toward easily accessible tissues like blood. For example, the eQTLGen consortium database provides largescale eQTL summary data from 31000 whole-blood samples, and analogous large database exist for blood-based mQTL.²³ Because of the lack of spine-specific QTL data, our study is only for exploratory in nature. We then conducted thorough functional enrichment analyses that included Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and Reactome pathways on the selected candidate genes. Furthermore, we conducted computational prediction of therapeutic small-molecule compounds through the Connectivity Map database followed by molecular docking simulations using AutoDock Vina to evaluate ligand-receptor binding affinities. This study aims to investigate how these genetic determinants contribute to the pathophysiology of spondylolisthesis development. Our goal is to identify molecular biomarkers for diagnosis and potential treatment of spondylolisthesis.

Materials and Methods

Summary of Databases, Tools, Chemicals and Instruments Used in the Study

A summary of the databases, tools, and other resources utilized in this study, along with the corresponding access links, is provided in Table 1.

Table 1 Summary of Databases, Tools, Chemicals and Instruments Used in the Study

Data Extraction

The dataset relevant to spondylolisthesis, designated as finn-b-M13_SPONDYLOLISTHESIS, was retrieved from the Integrative Epidemiology Unit Open genome-wide association study (IEU OpenGWAS) database. This dataset encompassed 16,380,280 single nucleotide polymorphisms (SNPs) derived from 167,351 European samples (cases: controls = 2,669: 164,682).

Sources of Exposure Factors

The expression quantitative trait Loci (eQTL) refer to genetic variations, often SNPs, that affect the expression levels of specific genes within an individual’s genomic landscape.²⁴ The study of eQTLs offers valuable insights into the genetic regulators of gene expression, shedding light on gene function and its association with phenotypic traits. The summary-level data are obtained from the eQTLGen consortium, which includes information on 10,317 SNPs associated with the traits of 31,684 individuals, primarily derived from blood samples of healthy Europeans. However, the eQTLGen consortium’s data exclude variants that are linked to the expression levels of genes situated on the X and Y chromosomes.²⁵ Furthermore, the associations between SNPs and cytosine-phosphate-guanine (CpG) sites were examined using blood DNA Methylation Quantitative Trait Loci (mQTL) data from 1,980 individuals of European ancestry, as reported by McRae et al.²⁶ The mQTL were predominantly composed of cis-mQTL.

Selection of Instrumental Variables (IVs)

The selection criteria for the IVs were as follows: 1) All SNPs incorporated within the primary analysis were required to have a p value of at least 5×10-8; 2) SNPs with an R² > 0.9 or R² < 0.05 in proximity to the leading SNPs were excluded, with only those pairs remaining where R² ≤ 0.9; 3) Linkage disequilibrium (LD) pruning in the GWAS data was carried out using information from the European population of the 1,000 Genomes Project. The heterogeneity in dependent instruments (HEIDI) approach was employed to ascertain whether the gene SNP-mediated phenotypes were influenced by LD response, with the criterion of HEIDI p value >0.05. The application of a HEIDI test p value below 0.05 indicated a heterogeneous association, suggesting the potential for pleiotropy. This step was implemented to guarantee the reliability and stability of the results.

Summary-Data-Based Mendelian Randomization (SMR) Analysis

SMR is a statistical technique that leverages genetic variation SNPs as IVs to infer causal relationships between exposures and outcomes, particularly in the study of complex diseases or traits, when direct randomized controlled trials are not feasible. The following methodology was employed to assess causality while ensuring adherence to MR analysis principles: βSMR = βSNP-GWAS/βSNP-QTL. In this formula, βSMR denoted the estimated effect size of eQTL/mQTL on spondylolisthesis GWAS, where βSNP-QTL was the estimated effect size of SNP on eQTL/mQTL (a genetic variant-exposure trait association) and βSNP-GWAS represented the estimated effect size of SNP on spondylolisthesis (the same genetic variant-outcome trait association). Which, a positive association was indicated by β > 0, while β < 0 suggested a negative association. The odds ratio (OR) was utilized to quantify the effect of an eQTL/mQTL on spondylolisthesis, calculated as OR = exp (βSMR). In this equation, OR represented the estimated ratio for each 1-ln increase at the genomic level of the eQTL/mQTL, with “exp” denoting the base of the natural logarithm. In this study, SNPs were utilized as IVs, with mQTLs and eQTLs serving as exposure factors and spondylolisthesis as the outcome. SMR analysis was conducted utilizing the SMR software (v 1.3.1) to identify genetic variants associated with spondylolisthesis. Specifically, multiple testing of the SMR results using the Bonferroni correction indicated that p values below 0.05 for eQTL or mQTL in relation to spondylolisthesis were considered statistically significant. Additionally, eQTLs or mQTLs with a p value <0.05 following SMR analysis and an HEIDI p value >0.05 were selected for further colocalisation analysis.

Population Scope and Gene Selection Criteria

In this study, analyses were confined to participants of European ancestry to maintain population homogeneity and minimize confounding from genetic stratification. The primary GWAS and QTL summary datasets available for spondylolisthesis were all derived from European cohorts. Focusing on a single ancestry ensures the validity of Mendelian randomization assumptions. Expansion to multi-ethnic data was not feasible due to lack of comparable datasets in other populations. Additionally, during candidate gene screening, we applied stringent filtering criteria to ensure biological relevance. Analyses were restricted to Homo sapiens genes, since all eQTL and mQTL data were human-derived, ensuring species-specific relevance.

Bayesian Colocalization Analysis

Colocalization analysis was utilized to genetically determine whether two potentially related phenotypes shared common genetic causal variants within a specified region. The Bayesian colocalization analysis was performed utilizing the coloc package (v 5.2.3) for the identified mQTLs or eQTLs (trait 1) and spondylolisthesis (trait 2).²⁷ The colocalization analysis was based on the following five exclusivity assumptions: 1) Hypothesis 0 (H0): no association with any traits; 2) H1: association only with trait 1; 3) H2: association only with trait 2; 4) H3: distinct causal variants for the two traits; and 5) H4: shared causal variants between the two traits. The posterior probabilities of H4 (PPH4) were examined by assessing the colocalization of all SNPs situated within a 100 kb region upstream and downstream of the lead SNPs, both in eQTL-GWAS and mQTL-GWAS analyses. A PPH4 value exceeding 0.7 was identified as an indicator of robust colocalization, and genes demonstrating this level of colocalization were designated as key genes due to their probable direct causal link with spondylolisthesis. After that, scatter plots were generated for the eQTLs and mQTLs of the key genes and spondylolisthesis using the SMR software (v 1.3.1). Of these, a positive slope in the scatter plot indicated that the eQTL or mQTL served as a protective factor against spondylolisthesis, whereas a negative slope suggested that it functioned as a risk factor for the condition.

Quantitative Reverse Transcription PCR (qRT-PCR) for Validation

Based on single-cell sequencing studies demonstrating that the human nucleus pulposus is predominantly composed of chondrolineage cells, such as regulatory chondrocytes, fibro-chondrocytes and pre-chondrocytes, the SW1353 chondrosarcoma cell line was chosen for this investigation.^28,29 SW1353 cells (Procell; China, Wuhan), a well-characterized human chondrosarcoma cell line, were employed as the experimental model. The treatment group received 500 μM hydrogen peroxide (H₂O₂) for 24 hours, while the control group was administered an equivalent volume of physiological saline under identical conditions.

RNA extraction was performed using a phenol-chloroform method, followed by quantification via ultraviolet spectrophotometry. Subsequent reverse transcription was conducted using the RevertAid First Strand cDNA Synthesis Kit (AE301-02 TransGen). qRT-PCR amplification was carried out with SYBR Green chemistry on a QuantStudio™ 6 Flex system (Applied Biosystems).

Gene expression levels were calculated using the 2^−ΔΔCT method, with normalisation to the endogenous reference gene GAPDH. Differential gene expression between groups was assessed using Student’s t-test, where P < 0.05 was considered statistically significant.

LOXL-4

Forward primer Sequence (5’->3’): ATGAGACCTGCCTCTTCCAC

Reverse primer Sequence (5’->3’): AGTTGAGATCGGACGTGGC

GAPDH

Forward primer Sequence (5’->3’): CTTTGCAGCAATGCCTCCTG

Reverse primer Sequence (5’->3’): GAAAGGTGGGAGCCTCAGTC

Tissue-Specific Expression and Subcellular Localization

The Human Protein Atlas (HPA) database was employed to investigate the expression patterns of crucial genes across different human tissues and regions of the brain, offering valuable information regarding the mRNA and protein levels of these genes. This study aided in examining potential mechanisms through which these key genes might act as therapeutic targets for spondylolisthesis treatment. Subcellular localization refers to determining the exact location within a cell where the product of a specific gene (usually a protein) is expressed, shedding light on the functional role of that product. The subcellular localization of key genes was assessed using the mRNALocater database.

RNA Methylation Modification and Protein Phosphorylation

The methylation modification sites of biomarkers were predicted to gain further insights into the regulatory mechanisms of key gene expression. The sequence-based RNA adenosine methylation site predictor (SRAMP) was utilised to predict N6-methyladenosine (m6A) modification sites on each pivotal gene, in conjunction with their specific locations within the RNA secondary structure. Furthermore, phosphorylation represents the most prevalent form of covalent modification within post-translational protein processes. The potential phosphorylation sites for each key gene were analyzed individually using the PhosphoSitePlus online tool.

Gene Ontology (GO) Enrichment Analysis

The biological functions of the key genes were elucidated through the use of the clusterProfiler package (v 4.7.1.003), which enabled the execution of GO enrichment analysis on the key genes (adjusted p value < 0.05).³⁰ Specifically, the GO comprised 3 principal sections: biological process (BP), cellular component (CC), and molecular function (MF).

lncRNA-Key miRNA-mRNA and TF-mRNA Interaction Network Building

It has been shown that molecular regulatory networks provide important understanding of the core processes involved in gene regulation during disease progression. To investigate the regulation of key genes by microRNAs (miRNAs), the TargetScan database was utilized in conjunction with the ElMMo and PITA databases for predicting miRNAs targeting these key genes. Subsequently, the intersection of the predicted miRNAs from the aforementioned three databases was determined using the VennDiagram package (v 1.7.3), yielding the key miRNAs.³¹ After that, the upstream long non-coding RNAs (lncRNAs) of key miRNAs were predicted utilizing the StarBase database, employing a screening criterion of clipExpNum >4. Moreover, the transcription factors (TFs) interacting with key genes were predicted utilizing the miRNet database. The regulatory network maps for lncRNA-key miRNA-mRNA and TF-mRNA were visualized using the Cytoscape software (v 3.5.0).³²

Prediction of Small Molecule Compounds and Molecular Docking

Prediction of small molecule compounds and molecular docking through the comparative toxicogenomics database (CTD), we identified small molecule compounds that may effectively target significant genes. A network map illustrating the interactions between these small molecule compounds and key genes was generated using Cytoscape software (version 3.5.0). Following this, we selected the three small molecule compounds that demonstrated the highest interaction counts with the target genes for subsequent molecular docking studies. To create the 3D structures of these key genes, we employed AlphaFold v2.0, an artificial intelligence-based prediction tool, for in silico modeling. The amino acid sequences of the relevant key genes were retrieved from the UniProt-KB database. We evaluated the stereochemical quality of these predictive models using the Local Distance Difference Test (LDDT) score provided by the AlphaFold database, where scores exceeding 90 were considered indicative of outstanding stereochemical accuracy for the key genes. We also used a gradient of green to represent the Predicted Aligned Error(PAE). Lighter green denotes low PAE values and high reliability in rigid regions, and darker green signifies high PAE values and uncertainty in flexible or poorly defined areas. Following protein structure analysis, molecular docking was utilized to simulate the optimal orientation of small molecule compounds binding to macromolecules, thereby facilitating a more profound comprehension of small molecule compounds-gene interactions through computational simulation. The candidate small molecule compounds were imported into the public chemistry database (PubChem) to obtain their 3D structures. Following this, the key genes were subjected to molecular docking with the candidate small molecule compounds utilizing the CB-Dock2 online tools. Since LOXL4 lacked a well characterized binding pocket, a blind docking method was employed. The CB-Dock2 platform was used to automatically identify the top-ranking potential cavities on the protein surface. The initial grip box dimensions were set to approximately 20Å per side and were appropriately adjusted according to the volume of each cavity. The docking procedure utilized the default parameters of Vina, with the global search exhaustiveness set to 8 and a maximum of 10 binding poses output per ligand. The results are presented as binding free energy estimates. Lower binding free energy values indicate a more stable ligand-receptor interaction. Based on previous studies, a blinding free energy ≤-5.0kcal/mol is generally considered indicative of a strong binding affinity between the ligand and the receptor.^33,34 Therefore, we adopted −5.0kcal/mol as the cutoff for this study.

Statistical Analysis

The statistical analyses were conducted utilizing R language (v 4.2.2), and the P value <0.05 was deemed to be statistically significant.

Results

eQTL, mQTL and SMR Analysis of Spondylolisthesis Outcomes

After HEIDI test (p value > 0.05) and Bonferroni correction (p value < 0.05), 841 cis-eQTLs and 3,224 mQTLs were identified, which were found to be significantly associated with spondylolisthesis (Figure 1A and B, Supplementary Tables 1 and 2). Subsequently, the Bayesian colocalization analysis revealed that the PPH4 values for both eQTL (ENSG00000138131) and mQTL (cg09335911) of LOXL4 surpassed 0.90 (Table 2). The results indicated a robust potential correlation between LOXL4 and spondylolisthesis, presumably attributable to shared genetic variants. Confirmation of the credibility of the hypothesis was provided by a higher posterior probability, which identified LOXL4 as a key gene influencing the progression of spondylolisthesis. Furthermore, the scatter plots for ENSG00000138131 and cg09335911, both positioned within the LOXL4, exhibited negative slopes, indicating that LOXL4 may serve as a risk factor for spondylolisthesis (Figure 1C and D).

Table 2 LOXL4 eQTL and mQTL Bayesian Co-Localization Analysis

Figure 1 The eQTL, mQTL and SMR analysis of spondylolisthesis. (A) Manhattan plot of eQTL vs SMR analysis of SPONDYLOLISTHESIS GWAS data with chromosome positions in horizontal coordinates and -log10(P.Value) in vertical coordinates, dashed line represents the threshold for significance after applying the Bonferroni correction and HEIDI test, indicates p-value of 0.05. (B) Manhattan plot of mQTL vs SMR analysis of SPONDYLOLISTHESIS GWAS data. Manhattan plot. Dashed line indicates p-value of 0.05. (C) Scatterplot of eQTL (SMR analysis of SPONDYLOLISTHESIS GWAS data, with the effect of eQTL or mQTL). The horizontal coordinate represents the effect size of SNP on eQTL, while the vertical coordinate indicates the effect on SPONDYLOLISTHESIS outcome. Dashed line indicates the fitted results derived from the SMR analysis. (D) Scatterplot of mQTL (SMR analysis of SPONDYLOLISTHESIS GWAS data, with the effect of eQTL or mQTL). The horizontal coordinate represents the effect size of SNP on mQTL, while the vertical coordinate indicates the effect on SPONDYLOLISTHESIS outcome. Dashed line indicates the fitted results derived from the SMR analysis.

LOXL4 Was Validated by qRT-PCR

SW1353 cells were utilised to establish experimental and control groups for investigating differential LOXL4 gene expression. qRT-PCR analysis demonstrated that hydrogen peroxide (H2O2) induced chondrocyte degeneration significantly upregulated LOXL4 expression (P < 0.01) (Figure 2). These findings are consistent with the hypothesis that LOXL4 responds to degenerative stimuli.

Figure 2 qRT-PCR Analysis of LOXL4 Expression upregulation in H2O2-Induced Degenerated human osteosarcoma cells (SW-1353).

Tissue-Specific Expression Patterns and Subcellular Localization of LOXL4

The distribution of LOXL4 in tissues was further investigated. As indicated by predictions derived from the HPA database, LOXL4 demonstrated the highest expression levels within the salivary gland and the choroid plexus of the brain (Figure 3A and B). This finding might provide novel insights into the target tissue and pathway of LOXL4. In addition, LOXL4 was predominantly localized in the cytoplasm (Figure 3C), implying a crucial role in its cytoplasmic synthesis, processing, and secretion.

Figure 3 Distribution of key genes in the human body. LOXL4 gene expression (vertical coordinate) across various human tissues and specific brain regions (horizontal coordinate). (A) Expression of key genes in human tissues. (B) Expression of key genes in the brain. (C) Subcellular localization map of key genes.

RNA Methylation Modification and Protein Phosphorylation of LOXL4

The study of mRNA methylation modification was crucial for elucidating gene expression regulation mechanisms. The results indicated that LOXL4 may possess numerous m6A methylation sites, with several of these sites being highly confident m6A methylation locations (Figure 4A). m6A is a prevalent methylation modification in mRNAs, and the identification of m6A methylation sites in LOXL4 could be significant for investigating its gene expression regulation mechanism. Furthermore, predictions from the PhosphoSitePlus database suggested that LOXL4 might be translated into two proteins, SRCR and Lysyl_oxidase, with the SRCR protein primarily undergoing phosphorylation modification (Figure 4B).

Figure 4 Modification site prediction. (A) Prediction of m6A sites of key genes, horizontal coordinate is the position, vertical coordinate is the modification score, and different color lines are the confidence level, red, purple, blue, and green are very high, high, medium, and low confidence level, respectively. (B) Post-translational modification sites of key genes, horizontal coordinate is the position, vertical coordinate is the number of times they have been reported in the literature.

Functional Analyses Related to LOXL4

The enrichment analysis was conducted to provide preliminary insights into the signaling pathways implicated by LOXL4. A total of 10 GO entries were found to be significantly enriched with the protein LOXL4, including six BPs, one CC, and three MFs. In particular, the analysis of pathways revealed the critical involvement of LOXL4 in the processes of “protein oxidation”, “collagen-containing extracellular matrix”, and “oxidoreductase activity” (Figure 5). These findings suggested that LOXL4 played a significant role in regulating the stability of the extracellular matrix, the structural repair of cells and tissues, and its antioxidant function.

Figure 5 GO enrichment analysis of key genes, horizontal coordinate is the number of genes, vertical coordinate is different pathways.

LOXL4 Was Regulated by Multiple Factors

The analysis of molecular regulatory networks provided a more comprehensive understanding of the regulatory factors that influence LOXL4. A total of 8 key microRNAs were identified by taking the intersection of the TargetScan, ElMMo and PITA datasets (Supplementary Figure 1). Subsequently, 53 lncRNAs were predicted to target the key miRNAs. Consequently, a lncRNA-key miRNA-mRNA network encompassing LOXL4, along with 8 key miRNAs and 53 lncRNAs, was established (Figure 6A), suggesting multi-factorial regulation of LOXL4. To illustrate, SNHG4 was observed to modulate LOXL4 expression through the regulation of hsa-let-7g-5p. Furthermore, LOXL4 was predicted to be regulated by 9 TFs, such as SOX5 (Figure 6B).

Figure 6 ceRNA network analysis. (A) mRNA-miRNA-lncRNA network diagram comprise 62 nodes with 431 edges, mRNAs in blue, miRNAs in purple, lncRNAs in pink. (B) mRNA-TF network diagram comprise 10 nodes with 9 edges, mRNAs in blue, TFs in pink.

Small Molecule Compounds Prediction and Molecular Docking of LOXL4

The interactions between small molecule compounds and LOXL4 were obtained from the CTD database, which resulted in the identification of 93 potential therapeutic agents for the treatment of spondylolisthesis (Figure 7A). After that, the crystal structure of LOXL4 was predicted utilizing AlphaFold v2.0, the LDDT score and low PAE indicated favorable stereochemical properties of the predicted structure (Figures 7B and C). Following an initial virtual screening based on interaction count, the top candidate compounds underwent further evaluation. From this pool, bisphenol A was excluded due to toxicity, the third and fourth compounds failed preliminary docking, and the fifth compound was omitted due to an unavailable structural formula. Consequently, estradiol and progesterone (the second and sixth compounds) were selected for final molecular docking analysis with LOXL4. The results showed strong predicted blinding affinities, with binding free energies of −8.4 kcal/mol for estradiol and −8.1 kcal/mol for progesterone, using a threshold of −5 kcal/mol to indicative high affinity (Figures 7D–G). These in silico findings suggest that sex hormones may interact directly with LOXL4, though further experimental validations is required to confirm this potential mechanism.

Figure 7 Drug prediction. (A) Key gene prediction drug map comprise 94 nodes with 93 edges, blue is gene, pink is drug. Each label in the map denotes a unique predicted drug, and all edges represent an association with the LOXL4. (B) 3D models of key genes, color is expected distance deviation, ranging from 0–31, darker color indicates better prediction. (C) Both horizontal and vertical coordinate represent the PAE values. Lighter green denotes low PAE values and high reliability in rigid regions, and darker green signifies high PAE values and uncertainty in flexible or poorly defined areas. (D) Molecular docking diagram of LOXL4 with Bisphenol A. Yellow dashed lines show hydrogen bonding forces, color is small molecule compounds, blue is protein structure, and black text shows amino acid residues and positions of small molecule compounds and proteins connected by hydrogen bonding. (E) Overall view of molecular docking of LOXL4 with Progesterone. (F) Enlargement of the molecular docking between LOXL4 and Estradiol. (G) Localized magnification of molecular docking of LOXL4 with Progesterone.

Discussion

Spondylolisthesis, a degenerative spinal disorder characterized by vertebral displacement leading to spinal instability and neural compression, demonstrates inheritable genetic susceptibility.^35,36 The study identified LOXL4 as a new risk gene for spondylolisthesis by integrated multi-omics analysis using an SMR method based on pooled data (P < 0.05), Bayesian colocalization confirmation (PPH4 > 0.90), and excluding the effect of multiplicity of effects (HEIDI P > 0.05), and confirmed for the first time that dysregulation of LOXL4 expression is genetically linked to the pathogenesis of lumbar spondylolisthesis. Spatial expression profiling via the Human Protein Atlas demonstrated predominant LOXL4 expression in salivary gland and cerebral choroid plexus tissues. Functional network analysis revealed its involvement in lncRNA-miRNA regulatory axes and transcription factor interactions, suggesting epigenetic mediation in disease progression. Utilizing the Comparative Toxicogenomics Database (CTD), we prioritized estradiol and progesterone as high-potency therapeutic agents targeting LOXL4, with molecular docking simulations confirming strong binding affinities (binding energies < −5.0 kcal/mol) through stable hydrogen bonding and hydrophobic interactions. These findings provide mechanistic insights into LOXL4-mediated spinal degeneration and propose repurposable endocrine modulators for pharmacological targeting.

Lysyl oxidase (LOX) and its related enzymes (LOXL1-4) play critical roles in extracellular matrix (ECM) homeostasis and tissue mechanical properties.³⁷ These enzymes catalyze the oxidative deamination of lysine residues, promoting covalent cross-linking between collagen and elastin, thereby enhancing the mechanical strength and elasticity of the ECM.³⁸ Within spinal architecture, type I collagen (predominant in the annulus fibrosus) provides tensile strength, while type II collagen (in the cartilaginous endplate) maintains hydrophilic properties for structural integrity. Simultaneously, elastin governs the elastic recoil of ligamentum flavum, ensuring dynamic stability during spinal movement and loading. In degenerative processes such as spondylolisthesis, abnormal collagen metabolism manifests as disc dehydration, annular fissures, and ligamentous laxity. These alterations may compromise spine stability and facilitate vertebral displacement.^38,39

From another perspective, the LOX/LOXL family not only contributes to structural support but may also modulate ECM cell interactions, influence matrix remodelling dynamics, affect cellular morphology, migration and apoptosis through variations in expression and enzymatic activity.⁴⁰ In the context of spinal degeneration, dysregulated LOXL4 expression could impair the required cross-linking density and stiffness of the local ECM, rendering it more susceptible to structural deformation or failure under mechanical load and fluctuating stress conditions.

As a member of the LOX family, LOXL4 has not been extensively studied in the context of spinal disorders. However, its potential role in collagen cross-linking and tissue homeostasis warrants attention. In a model of pulmonary fibrosis, LOXL4 has been identified as a key mediator of pathological collagen cross-linking and proposed as a potential therapeutic target.⁴¹ Furthermore, in other research, a family-based whole-exome sequencing study suggested that LOXL4 variants may influence bone-cartilage mechanical properties by altering collagen cross-linking and bone matrix structure. (Family-based whole-exome sequencing implicates a variant in lysyl oxidase like 4 in atypical femur fractures).

Notably, our qRT-PCR validation experiments showed that oxidative stress significantly upregulated LOXL4 in SW1353 cells, consistent with the idea that LOXL4 is responsive to degenerative stimuli. However, these results are based on a single in vitro chondrosarcoma cell model and cannot fully capture the complex in vivo context of spondylolisthesis. In particular, the facet joints of the spine, which are composed of cartilage, play an important role in spinal stability and in the pathogenesis of spondylolisthesis. Future studies should therefore investigate LOXL4 expression and function in facet joint cartilage under degenerative conditions.

As a vital stabilizing component of the spine, degeneration of the intervertebral disc is closely linked to spondylolisthesis. The pathological alterations associated with intervertebral disc degeneration (IVDD) encompass oxidative stress, apoptosis, and inflammatory factors.⁴² In the functional enrichment analysis conducted in this study, it was found that LOXL4 plays a significant role in pathways including “protein oxidation,” “oxidoreductase activity,” and “collagen-containing extracellular matrix (ECM).” Studies have demonstrated that oxidative stress can trigger intervertebral disc degeneration through multiple pathways; consequently, strategies that inhibit the overproduction of reactive oxygen species (ROS) within the disc while enhancing their removal have been shown to effectively postpone disc degeneration.^43–46 This process may also be related to the pathological mechanisms involved in spondylolisthesis.

Oxidoreductases are also involved in ROS metabolism. The delicate balance between ROS and antioxidants is crucial for maintaining the normal function and survival of intervertebral disc cells. Excessive levels of ROS levels can damage macromolecules such as nucleic acids, lipids, and proteins, affecting normal cellular activities and functions, ultimately leading to cell aging or death.^47–49 Recent evidence demonstrates that miR-193a-3p directly targets LOXL4 to modulate oxidative stress pathway activity, thereby potentiating multidrug chemoresistance in bladder carcinoma cells.⁵⁰ This finding suggests LOXL4 may mediate diverse pathological processes through regulation of cellular redox homeostasis. Similarly, under pathological mechanical loading (eg, muscular traction), spinal stabilization structures such as paraspinal muscles and facet joints exhibit oxidative stress-mediated activation of TGF-β signaling pathways. These molecular events degrade the fibroelastic integrity of the tissues, ultimately elevating risks of vertebral displacement.⁵¹

The structural integrity of intervertebral discs and ligaments is critically dependent on ECM homeostasis, which is maintained through coordinated interactions among elastic fibers, collagen, and aggrecan.^52,53 However, disc degeneration involves disruption of anabolic-catabolic equilibrium.^54,55 LOXL4 enzymatically catalyzes collagen-elastin crosslinking, directly modulating ECM structural stability.^56,57 Histopathological analyses reveal three characteristic ECM abnormalities in degenerative spondylolisthesis: elastic fiber depletion, collagen network fragmentation, and pathological ECM remodeling.^58,59 These findings implicate LOXL4 dysfunction in destabilizing intervertebral joints through impaired interfibrillar cross-linking. Emerging evidence further indicates that dysregulated LOXL4 expression disrupts ECM remodeling cascades, contributing to multisystem pathophysiology.⁵⁶ Mechanistically, LOXL4 drives fibrotic progression via dual mechanisms: direct ECM restructuring through pathological collagen crosslinking, and activation of pro-fibrotic signaling pathways (eg, YAP/TAZ and TGF-β). This establishes a self-perpetuating fibrotic microenvironment.⁶⁰ Analogous molecular mechanisms underlie hypertrophic remodeling of spinal ligaments.^61,62

Molecular docking analysis reveals high-affinity binding of estradiol and progesterone to LOXL4 (binding energy < −5.0 kcal/mol), suggesting that direct ligand-enzyme interactions may modulate LOXL4 bioactivity. However, this remains a hypothesis. Thus, while sex hormones may influence LOXL4, it would be premature to propose their use as therapeutic modulators without further evidence. Intriguingly, beyond classical reproductive functions, these steroid hormones regulate connective tissue homeostasis through genomic (nuclear ER/PR) and non-genomic (membrane receptor) pathways.^63–65 In models of intervertebral disc degeneration, estradiol suppresses collagenolysis via ERα-mediated MMP-13 downregulation, whereas progesterone enhances TIMP-1 production through PR-B signaling to preserve ECM equilibrium.^66–68 LOXL4 enzymatic activity may be coregulated by hormone-receptor complex formation, suggesting crosstalk between endocrine signaling and ECM dynamics.

Degenerative spondylolisthesis demonstrates significant sexual dimorphism in adults, which is potentially mediated by sex hormone fluctuations.^69,70 A systematic review by Marcin et al⁷¹ established causal links between sex hormone deficiency and musculoskeletal remodeling, which significantly elevates spondylolisthesis risk. However, epidemiological discrepancies persist regarding the estrogen-spondylolisthesis association.⁷² Cross-sectional data indicate oophorectomy-induced menopause significantly increases sarcopenia prevalence.⁷³ A longitudinal follow-up with an average duration of 19.8 years indicated that estrogen replacement therapy in women who are postmenopausal improves vertebral bone mineral density (BMD) and decreases the occurrence of spinal pathologies, implying that there are osteoprotective mechanisms linked to the stabilization of the musculoskeletal system.⁷⁴ Conversely, a 4-year longitudinal study by Wáng et al⁷⁵ in Chinese elderly populations found no direct correlation between estrogen levels and spondylolisthesis incidence. Using ovariectomized murine models, Hao et al⁷⁶ documented accelerated zygapophyseal joint cartilage degeneration under estrogen-deficient conditions. Pathological reduction of estrogen receptor alpha (ERα) in articular chondrocytes induces osteoarthritic changes.⁷⁷ These endocrine-cartilage interactions may constitute key etiological factors for postmenopausal spondylolisthesis.⁷⁸

Our findings suggest a potential connection between LOXL4 and sex hormones. However, it is important to note that the GWAS, eQTL, mQTL datasets utilized in this study were not sex-balanced, nor did they include any sex-stratified analyses. The potential for sex-biased genetic effects to go undetected in combined sex analyses has been noted in few GWAS literatures, meaning potential sex-specific genetic effects are often not explored.⁷⁹ But there is evidence that biological sex can influence gene regulation, for example, estradiol has been shown to upregulate the expression of LOX family genes in certain tissues.⁸⁰ Some genetic associations are known to differ by sex.⁸¹ In this study, the absence of explicit sex-specific date in our analysis limits the interpretation of the link between LOXL4 and sex hormone. The observation of this research is exploratory and requires additional studies. Specifically, future research should validate whether modulating LOXL4 with sex hormones has any therapeutic significance in spondylolisthesis.

Limitation

This study has some limitations. Firstly, the gene expression and methylation data utilized in this study were derived primarily from blood samples in publicly available eQTL and mQTL datasets. While these datasets provide a broad overview of genetic associations. These may not fully reflect the tissue-specific expression patterns of LOXL4 in the spine and the pathology of spondylolisthesis manifests. This discrepancy could introduce potential bias, as gene regulation mechanisms in the intervertebral disc and facet joints may differ from those observed in blood. In the absence of spine-specific QTL data, any causal inferences drawn must be interpreted with caution given the potential limitations imposed by tissue-specific expression differences. Secondly, although the HEIDI test was employed to mitigate the influence of LD, it does not fully eliminate the risk of horizontal pleiotropy. Future studies that incorporate advanced techniques to detect and correct for pleiotropy will be important to validate our findings and avoid potential confounding. Thirdly, our in vitro validation was limited to a single cell line and a single outcome measure. This study examined LOXL4 expression changes in the SW1353 cell line under oxidative stress, but this immortalized cell model may not fully recapitulate the behavior of primary nucleus pulposus or annulus fibrosus cells in an actual degenerative disc environment. Moreover, focusing solely on LOXL4 mRNA as the endpoint provides a narrow view of the degenerative process. We did not assess downstream effects such as collagen degradation, cell viability and proteoglycan loss. These factors limit the generalizability of our experimental findings. Lastly, while the interaction between LOXL4 and sex hormones was predicted. It remains unclear whether this interaction plays a direct role in spondylolisthesis progression.

Conclusion

In conclusion, this exploratory study employed an integrative approach that combined SMR, bayesian colocalization, molecular docking, and in vitro qRT-PCR validation to investigate the potential association of LOXL4 in degenerative spondylolisthesis. Our analyses yield preliminary, multi-faceted evidence suggesting an association between LOXL4 dysregulation and disease pathogenesis, which may involve its known functions in extracellular matrix biology and oxidative stress response. Notably, molecular docking simulations propose a speculative yet intriguing interaction between LOXL4 and sex hormones (estradiol and progesterone), framing a novel hypothesis for hormonal influence on ECM homeostasis. It is crucial to emphasize that the hormonal link remains speculative, and the functional role of LOXL4 inferred from our in vitro model requires validation in more physiologically relevant systems. Therefore, the precise functional contribution of LOXL4 to spondylolisthesis in vivo, along with its proposed endocrine modulation, remains to be established. Collectively, this work outlines a plausible mechanistic framework and highlights LOXL4 as a candidate warranting further investigation. Future studies should focus on validating LOXL4 function in spinal tissues and elucidating whether its activity can be modulated therapeutically, including testing the proposed interaction with sex hormones.