Jianguo Zhou1†
Jian Kang- 1Orthopedics and Traumatology, The Second Hospital of Liaoning University of Traditional Chinese Medicine, Shenyang, Liaoning, China
- 2Orthopedics and Traumatology, Affiliated Hospital of Liaoning University of Traditional Chinese Medicine, Shenyang, Liaoning, China
- 3Phase I Clinical Research Unit, Affiliated Hospital of Liaoning University of Traditional Chinese Medicine, Shenyang, Liaoning, China
- 4First Clinical College, Liaoning University of Traditional Chinese Medicine, Shenyang, Liaoning, China
Background: While the association between hip osteoarthritis (HOA) and femoral neck bone mineral density (FN-BMD) is established, their shared genetic architecture remains elusive. This study aims to explore the genetic correlation and underlying mechanisms.
Methods: The present study applied bidirectional Mendelian randomisation (MR) to investigate causal relationships between HOA and FN-BMD. The quantification of genetic correlations was achieved by employing linkage disequilibrium score regression (LDSC) and high-definition likelihood (HDL) inference. Multi-trait analysis of genome-wide association studies (MTAG) have been shown to enhance statistical resolution, thereby unveiling hitherto unreported genetic associations. Independent MTAG-HOA risk loci were identified through conditional joint analysis (GCTA-COJO), complemented by functional mapping and annotation (FUMA) functional annotation. The application of both MAGMA and GCTA-fastBAT has revealed pleiotropic genes linked to MTAG-HOA susceptibility. Integration of fine-mapped genes from these loci with risk-associated candidates has enabled the identification of 13 key HOA-related genes. Functional annotation of these 13 key genes was performed using Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses. Multi-tissue transcriptome-wide association studies (TWAS) explored the expression of key genes across different tissues and their association with HOA. SMR analysis evaluated the causal relationship between key gene expressions in various tissues and HOA. Proteomic profiling is conducted via proteome-wide association studies (PWAS) and biomarker level imputation from summary statistics (BLISS). The application of stratified LDSC-SEG has revealed a genetic enrichment profile in cell types.
Results: Bidirectional MR analysis revealed a significant negative causal effect of FN-BMD on HOA (β = −2.17, P < 0.01), whereas the reverse MR analysis did not identify a causal effect. LDSC and HDL analyses revealed genetic correlations between HOA and FN-BMD of rg = 0.132 and rg = 0.1697, respectively. GCTA-COJO and FUMA collectively identified 28 independent risk SNPs associated with HOA. MAGMA and GCTA-fastBAT identified 48 pleiotropic genes. Integrating independent risk loci and pleiotropic genes culminated in the identification of 13 key genes associated with HOA. An enrichment analysis revealed that 13 key genes were significantly associated with biological processes integral to cartilage development, osteogenesis, cell proliferation, apoptosis, and stem cell differentiation. Multi-tissue TWAS and SMR analyses indicated that seven genes were associated with HOA across 22 tissues, with brain tissues accounting for 28.6%. Furthermore, PWAS and BLISS methods were utilized to analyze the proteomic features of these key genes. LDSC-SEG analysis revealed enrichment of HOA heritability in Cartilage, Lymphocytes, Oocytes, B Lymphocytes, Germ Cells, Osteoblasts, and Embryoid Bodies.
Conclusion: This study provides a comprehensive analysis of the genetic correlation between HOA and FN-BMD, elucidating shared genetic architecture and pinpointing key genes. These findings offer novel insights into the interplay between HOA and FN-BMD and highlight potential therapeutic targets.
1 Introduction
Osteoarthritis (OA) is a joint disease characterized by articular cartilage degeneration, subchondral bone remodeling and joint space narrowing (Motta et al., 2023). Its pathogenesis results from the intricate interplay between systemic predisposing factors and local biomechanical influences. Established risk factors reported include age, obesity, mechanical stress, metabolic disorders and genetic predisposition (Carpenter et al., 2006). Hip and knee joints, the primary weight-bearing joints in humans, are most affected by OA (Husted, 2012; Cui et al., 2020). Historically, OA has been regarded primarily as a disease of cartilage, with cartilage degradation as the central pathological feature. However, recent research has increasingly recognized the critical role of subchondral bone in the pathogenesis of OA. In hip osteoarthritis (HOA), femoral bone abnormalities manifest as both focal changes, including marginal osteophytes and subchondral bone cysts, and spatially heterogeneous changes, such as subchondral plate thickening and reduced trabecular bone mineralization (Chiba et al., 2014; Wang et al., 2016). Current understanding of changes in femoral bone density and microstructure in HOA focuses primarily on biomechanical stress distribution and biological responses (Auger et al., 2022). Nevertheless, the underlying mechanisms, particularly the role of genetic factors, remain to be fully elucidated.
It is estimated that the heritability of OA exceeds 50%. Previous genetic investigations have identified 27 risk loci across hip, knee, and hand osteoarthritis (Zengini et al., 2018), but locus overlap among these sites is limited. Mendelian randomization (MR) studies have evidenced a causal link between HOA and femoral neck bone mineral density (FN-BMD) (Qu et al., 2023). However, the genetic mechanisms underlying FN-BMD variation in HOA patients, and the extent of shared genetic architecture between HOA and FN-BMD, remain poorly understood. Genome-wide association studies (GWAS) have emerged as a robust tool for dissecting genetic risk factors in complex diseases. To address these knowledge gaps, the present study employs a suite of advanced post-GWAS analytical methodologies to systematically investigate the shared genetic architecture between HOA and FN-BMD, and to further identify potential pleiotropic genes. The overarching objective of this study is to comprehensively dissect the genetic architecture of FN-BMD variation in HOA, thereby providing a foundation for a deeper understanding of HOA pathogenesis and the development of targeted therapeutic interventions.
Research that focuses on just one disease might not identify the important genetic loci and molecular regulatory mechanisms. Therefore, multi-trait analytical approaches are essential for broadening the phenotypic spectrum under investigation, identifying associated risk loci and elucidating the shared genetic aetiology of different diseases (Lee et al., 2021). Shared genetic aetiology suggests potential pleiotropy, which often represents genetic confounders that link these traits (Sivakumaran et al., 2011; Verbanck et al., 2018). Consequently, cross-trait analyses leveraging correlations between GWAS data for HOA and FN-BMD are recommended to explore pleiotropic genetic variants or loci across multiple traits. These pleiotropic loci could be potential therapeutic targets, providing opportunities for simultaneous prevention or treatment of these conditions.
This study aimed to elucidate the genetic basis of HOA and FN-BMD, and to further investigate the biological mechanisms underlying HOA. Initially, bidirectional MR analysis was employed to assess the causal relationship between HOA and FN-BMD; genetic correlations were quantified using linkage disequilibrium score regression (LDSC) and high-definition likelihood (HDL) inference, and multi-trait analysis of GWAS (MTAG) was utilized to enhance statistical power and uncover novel genetic associations. Subsequently, independent HOA risk loci were identified and fine-mapped using conditional and joint analysis (GCTA-COJO) and FUMA (Functional Mapping and Annotation) for functional annotation. MAGMA and GCTA-fastBAT were then applied to identify pleiotropic genes associated with HOA susceptibility. Through the integration of genes mapped from risk loci and pleiotropic genes, 13 key HOA-associated genes were identified. Subsequently, functional annotation of these 13 key genes was performed using Gene Ontology (GO) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses. Following this, a multi-tissue Transcriptome-Wide Association Study (TWAS) was conducted to analyze the expression patterns of these key genes across various tissues and to assess their association with HOA. Summary-data-based SMR analysis was then used to further investigate the causal effects of gene expression in different tissues on HOA. Furthermore, the proteomic profiles of the 13 key genes were comprehensively analyzed using Proteome-Wide Association Study (PWAS) and the Biomarker Level Imputation from Summary Statistics (BLISS) method. Finally, stratified LDSC-SEG (LD Score Regression for specifically expressed genes) analysis was applied to identify cell-type-specific genetic enrichment, offering novel insights into the mechanisms of HOA (Figure 1).
Figure 1. Flow chart of this study.
2 Materials and methods
2.1 Data source
Summary statistics for HOA were obtained from a large GWAS available in the GWAS Catalog (GCST007091). This GWAS integrated genetic data from 393,837 individuals of European ancestry, including 15,704 HOA cases and 378,169 controls. Case samples were derived from the UK Biobank and the Arthritis Research UK Osteoarthritis Genetics Consortium (arcOGEN), with diagnoses based on clinical criteria (joint replacement) or radiographic assessment (Kellgren-Lawrence grade ≥ 2) (Tachmazidou et al., 2019).
Genetic variants associated with BMD were derived from the summary statistics of the largest publicly available GWAS meta-analysis conducted by the Genetic Factors for Osteoporosis Consortium (GEFOS). This dataset includes GWAS results for FN-BMD in 49,988 European ancestry, with BMD measurements obtained via dual-energy X-ray absorptiometry (DXA). The original GWAS analysis assessed the additive effects of single nucleotide polymorphisms (SNPs) with a minor allele frequency (MAF) > 0.5% on FN-BMD, accounting for key covariates such as sex, age, and body mass index (BMI). To facilitate downstream analyses, GEFOS applied weighting and standardization to the BMD phenotype data (mean = 0, standard deviation = 1) (Zheng et al., 2015).
2.2 Bidirectional MR analysis
MR analysis was employed to investigate the causal relationship between HOA and FN-BMD using SNPs as instrumental variables (IVs). The selection of IVs was based on fulfilling three core assumptions: 1 Relevance assumption: IVs must be strongly associated with the exposure; 2 Independence assumption: IVs must be independent of confounders; 3 Exclusion restriction assumption: IVs can only affect the outcome through the exposure and not through other pathways (Palmer et al., 2012).
The screening criteria for IVs included a strong association with exposure (p < 5 × 10−8), low linkage disequilibrium (LD) (r2 < 0.001), and F-values greater than 10 to mitigate weak instrument bias. Steiger filtering was applied to ensure the correct directionality of the causal effect. Outlier IVs were identified and excluded using the MR-PRESSO method. Potential confounders, such as estrogen, were identified and removed using the LDlink database (Lin et al., 2020). The primary MR analysis method used was Inverse Variance Weighted (IVW) analysis. Heterogeneity among IVs was assessed using Q-tests, and pleiotropy was evaluated using MR-Egger regression based on its intercept estimates.
2.3 Genetic evidence analysis
2.3.1 Linkage disequilibrium score regression (LDSC)
LDSC is a statistical method used to estimate the genetic correlation (rg) between HOA and FN-BMD. LDSC evaluates the strength of the association between each SNP and the traits by estimating the LD score for each SNP. The results of cross-trait LDSC analyses are reported as genetic correlations and their standard errors (se). However, LDSC analysis may not provide reliable results if the heritability (h2) of either or both traits is low, which could limit the interpretation of the genetic correlation estimates (Bulik-Sullivan et al., 2015).
2.3.2 High-definition likelihood inference (HDL)
HDL regression was employed to enhance the precision and robustness of genetic correlation estimates and to validate the findings from LDSC analysis. Compared to LDSC, HDL offers a more advanced approach to genetic correlation analysis, significantly improving estimation accuracy. By leveraging leading eigenvalues with greater explanatory power and eigenvectors of the LD score matrix, HDL more effectively captures genetic correlation signals. HDL regression reduces the variance in genetic correlation estimates by approximately 60%, thereby increasing statistical power comparable to a substantial expansion in sample size (Ning et al., 2020). The reference dataset for HDL analysis included an LD matrix and its eigenvalue decomposition derived from 335,265 UK Biobank participants, incorporating 1,029,876 high-quality HapMap3 SNPs.
2.4 Multi-trait analysis
MTAG is a method that employs generalized inverse-variance-weighted meta-analysis to integrate summary statistics from multiple GWAS datasets. This method has been extensively applied to investigate the shared genetic architecture among traits, and it is currently recognised as a principal approach for exploring the common genetic underpinnings of diverse phenotypes (Guo et al., 2022; Liu et al., 2024; Meng et al., 2025; Song et al., 2025). The primary objective of MTAG is to enhance the statistical power of GWAS and improve the detection of associated loci by capitalizing on genetic correlations among traits. While the primary objective of MTAG-derived GWAS results is to strengthen signals for the primary trait, loci associated with this primary trait that exhibits collinearity with other traits will also be evident. This will aid in the discovery of potential shared (pleiotropic) loci (Turley et al., 2018). Furthermore, as this study aims to elucidate the genetic architecture of HOA, GWAS data for HOA and FN-BMD were combined to construct MTAG datasets (MTAG-HOA and MTAG-FN-BMD), and subsequently, the MTAG-HOA dataset was selected for further analyses.
2.5 Identification of independent risk loci
2.5.1 FUMA genetic association localization and annotation
The FUMA platform used stringent criteria with a maximum SNP P of < 5 × 10−8 and an additional inclusive significance threshold of P < 0.05. Independent SNPs were identified using an r2 threshold of < 0.6, with lead SNPs restricted to an r2 of < 0.1 within a 1 Mb radius. Genomic risk loci were defined by merged regions where lead SNPs were within 250 kb (Watanabe et al., 2017).
2.5.2 Conditional and joint multiple-SNP analysis (COJO)
GCTA-COJO was conducted using GCTA software. COJO identified SNPs independently associated with MTAG-HOA susceptibility, reducing redundancy from multiple SNPs in high LD that might otherwise represent the same underlying signal. That is, any two loci with less than 1 Mb on the same chromosome were analyzed conditionally to observe whether other loci were also associated with the phenotype. If a SNP lost statistical significance after conditional analysis, its association with the phenotype was explained by other SNPs and was excluded from further analysis. The stepwise model selection framework implemented in GCTA-COJO was applied to perform conditional and joint association analysis (Yang et al., 2012). This approach identified significant, independently associated variants within established genomic risk loci (P.mtag < 5 × 10−8) and validated additional signals using a joint P threshold of < 5 × 10−8. This study used the 1000 Genomes Project Phase 3 European ancestry reference dataset (Auton et al., 2015).
2.6 Genetic insights into MTAG-HOA
2.6.1 MAGMA analysis
MAGMA employs a multiple regression model to compute the cumulative effect of SNPs assigned to each gene (±10 kb). Subsequently, to quantify the gene-phenotype association, MAGMA aggregates SNP-level association statistics to derive gene scores (de Leeuw et al., 2015). The genome-wide Bonferroni-corrected significance level was set as P < 2.77 × 10−6 (0.05/18,079).
2.6.2 GCTA-FastBAT analysis
GCTA-fastBAT framework implements efficient gene-based association testing through quadratic form aggregation of SNP-derived z-statistics within defined genomic regions. This set-based approach calculates P using the approximate distribution of χ2 statistic summation, with LD correction via the 1,000 Genomes Project Phase 3 reference panel (Bakshi et al., 2016). The genome-wide Bonferroni-corrected significance level was set as P.FastBAT < 2.05 × 10−6 (0.05/24,390).
2.7 Functional enrichment analysis of MTAG-HOA-related key genes
Integration of findings from genes fine-mapped by GCTA-COJO and lead SNPs with those from pleiotropic genes identified through MAGMA and GCTA-fastBAT analyses consistently identified a specific set of genes as key determinants for MTAG-HOA. GO and KEGG databases were utilized for functional annotation, with a significance threshold of P < 0.05 applied. GO includes three domains: biological processes (BP), cellular components (CC), and molecular functions (MF) (Yu et al., 2012).
2.8 Multi-tissue transcriptome-wide association studie
The decision to incorporate all 49 tissues from GTEx v8 in the TWAS analysis was informed by the mounting recognition that skeletal health is not solely determined by local joint pathology, but also involves intricate interactions with the central nervous system (Zhang et al., 2025), peripheral nerves (Yoshino et al., 2021), digestive system (Liu et al., 2022), adipose tissue (Zapata-Linares et al., 2025), and immune responses (Xu et al., 2024). The GTEx v8 dataset provides comprehensive, high-quality expression reference panels across these major physiological systems, enabling tissue-specific evaluation of genetically regulated gene expressions (Li et al., 2018). Utilizing this comprehensive approach facilitates a more comprehensive investigation of the transcriptional architecture underlying MTAG-HOA and its genetic correlation with FN-BMD. By incorporating diverse tissue contexts, this study enhances the capacity to identify functionally relevant genes with context-dependent effects, thereby facilitating deeper insights into the multi-systemic etiology of hip osteoarthritis. MetaXcan leverages pre-established gene expression prediction models, which were trained within the GTEx v8 reference population and adjusted for sex, as well as potential technical and population-based confounding factors. These models facilitate the accurate estimation of the heritable component of gene expression levels in large, independent cohorts. The GTEx v8 prediction models are constructed upon the MASHR framework and, notably, do not simply incorporate all SNPs within a fixed window. Instead, they prioritize SNPs with a higher probability of causal involvement in gene expression regulation, using the DAP-G algorithm as a feature selection method (Wen et al., 2016). The DAP-G algorithm computes a posterior inclusion probability (PIP) for each SNP, retaining those with PIP >0.01 for model building (Urbut et al., 2019). This strategy is designed to improve both the predictive performance and biological interpretability of the models. As an extension of PrediXcan, MetaXcan integrates cis-eQTL weights derived from GTEx v8 and GWAS summary statistics, enabling analyses without requiring individual-level genotype data. This methodology presents a considerable advantage for the efficient identification of tissue-specific genes associated with traits of interest, particularly considering the growing availability of large-scale GWAS summary datasets. MetaXcan principally consists of two modules: S-PrediXcan and S-MultiXcan. S-PrediXcan was used to evaluate gene-trait associations within individual tissues, while S-MultiXcan was used to integrate association results across tissues (Barbeira et al., 2018; Barbeira et al., 2019). For all analyses, P-value were subjected to Benjamini-Hochberg (BH) correction, and a BH-adjusted P-value (P.adjust) < 0.05 was adopted as the threshold for statistical significance.
2.9 SMR analysis
SMR analysis leverages summary statistics from GWAS and eQTL to investigate potential pleiotropic associations between gene expression levels and complex traits. The Heterogeneity in Dependent Instruments (HEIDI) test is applied to evaluate potential horizontal pleiotropy within significant colocalization signals. The null hypothesis of the HEIDI test is the absence of horizontal pleiotropy within the colocalized signal. SMR and HEIDI methods together help determine whether alterations in gene expression levels mediate the SNP-phenotype effect, and to distinguish if this association reflects genuine causal mediation or is driven by linkage disequilibrium or horizontal pleiotropy (Zhu et al., 2016). This study employed SMR to investigate whether the expression levels of key genes mediate the genetic effects on MTAG-HOA, utilizing eQTL data from GTEx V8 as proxies for genetic variation. BH-adjusted P-value (P.adjust) < 0.05 was considered indicative of a significant SMR association. For the HEIDI test, P-value (P. HEIDI) > 0.01 suggests a lack of significant evidence for horizontal pleiotropy, thus supporting the interpretation that the SMR result potentially reflects a causally mediated effect of gene expression on the phenotype.
2.10 Proteomic studies
2.10.1 PWAS analysis
PWAS were performed using FUSION with publicly available, pre-calculated, protein-level prediction models from the INTERVAL study (Sun et al., 2018). These models were derived from the INTERVAL study, in which the levels of plasma proteins were measured in 3,301 individuals of European ancestry using the SomaLogic platform. Subsequently, 3,222 predictive models were developed for 3,170 proteins. The PWAS analysis made use of European linkage disequilibrium (LD) reference data from the 1000 Genomes Project.
2.10.2 BLISS analysis
BLISS is a novel method developed to create protein imputation models based on summary-level pQTL data, to conduct multi-ancestry PWAS analysis. By integrating proteomic data with GWAS findings, PWAS provides a more direct understanding of how genetic variations influence diseases through changes in protein abundance and function (Shan et al., 2025). In this study, drawing upon the deCODE study pQTL datasets, BLISS technology was employed to delineate the intricate proteomic signatures associated with MTAG-HOA.
2.11 LDSC-SEG
LDSC-SEG is an extension of stratified LD score regression (S-LDSC) designed to pinpoint relevant cell types or tissues for complex traits. The premise of LDSC-SEG is that if genomic regions associated with the heritability of a disease or trait are enriched for genes specifically expressed in a particular tissue, this tissue is likely to play a causal role in the disease etiology or trait biology. Utilizing gene expression profiles across 152 cell types curated by the Franke Lab as genomic annotations, LDSC-SEG assesses the expression levels of genes within these 152 cell types and ranks genes based on their cell-type specificity (Finucane et al., 2018).
3 Results
3.1 Bidirectional MR analysis
MR analysis, employing IVs selected based on three core assumptions, F-statistics, and the Steiger filtering method, revealed a negative causal effect of FN-BMD on HOA (β = −2.17,P < 0.01). The Egger intercept analysis detected no significant horizontal pleiotropy (P > 0.1), suggesting that the observed causal effect was not biased by such pleiotropy. Cochran’s Q test revealed no statistically significant heterogeneity across instrumental variables (P > 0.05), indicating consistency in the causal estimates derived from different instruments. Conversely, reverse MR analysis indicated no causal effect of HOA on FN-BMD (β = −0.091,P = 0.121). The Egger intercept analysis for the reverse direction also detected no significant horizontal pleiotropy (P > 0.1), suggesting the absence of bias from such pleiotropy in this analysis. In the reverse MR analysis, Cochran’s Q test indicated significant heterogeneity (P < 0.05), suggesting variability in causal estimates across the different instrumental variables (Figure 2).
Figure 2. Bidirectional MR forest map.
3.2 Genetic correlation between HOA and FN-BMD
Univariate LDSC analysis revealed that the h2 of HOA was 0.0325 (se = 0.0024), while that of FN-BMD was 0.2214 (se = 0.0249). Bivariate LDSC analysis further demonstrated a positive genetic correlation between HOA and FN-BMD (rg = 0.132, se = 0.0478, P = 0.0057).
HDL analysis revealed the heritability of HOA was 0.0307 (se = 0.0019), while FN-BMD exhibited the heritability of 0.1508 (se = 0.0169). In addition, there was a statistically significant positive genetic correlation between HOA and FN-BMD (rg = 0.1697, se = 0.0401, P = 2.36 × 10−5). Crucially, this positive genetic correlation robustly corroborates the findings from LDSC analysis, providing convergent evidence for a shared genetic architecture underlying HOA and FN-BMD (Table 1).
Table 1. Results of LDSC and HDL analyses.
Leveraging shared genetic signals between HOA and FN-BMD, MTAG boosts GWAS statistical power, resulting in an enhanced MTAG-HOA GWAS dataset. MTAG-HOA demonstrates a reduced false discovery rate (FDR) of 1.03%, indicating enhanced robustness of signal detection. Ultimately, the MTAG-HOA dataset comprises 6,673,816 validated SNPs, thereby conferring enhanced statistical power and reliability for subsequent analyses.
3.3 Risk loci for MTAG-HOA
The FUMA platform was utilized to identify lead SNPs for GWAS-FN-BMD, GWAS-HOA, MTAG-FN-BMD, and MTAG-HOA. Initial analysis identified 7, 44, 27, and 36 lead SNPs, respectively, representing independent genomic loci associated with GWAS-FN-BMD, GWAS-HOA, MTAG-FN-BMD, and MTAG-HOA. To further refine these signals and dissect independent genetic contributions, GCTA-COJO analysis was applied to the summary statistics of GWAS-FN-BMD, GWAS-HOA, MTAG-FN-BMD, and MTAG-HOA. This analysis identified 7, 14, 29, and 32 SNPs, respectively, which collectively accounted for independent association signals. Cross-validation of the FUMA-derived lead SNPs and the GCTA-COJO-selected independent SNPs yielded a robust set of 6, 6, 23, and 28 variants, respectively. The present study concentrated principally on the independent genetic risk loci for MTAG-HOA. The 28 SNPs were thus designated as independent genetic risk variants for MTAG-HOA, thereby further underscoring the complex polygenic architecture of HOA and highlighting key genetic loci contributing to disease susceptibility (Figure 3; Supplementary Tables S1, S2).
Figure 3. Venn diagram illustrating the intersection of GCTA-COJO and LeadSNP, and corresponding Manhattan plots with these risk loci annotated. (A,B) Venn diagram and Manhattan plot for FN-BMD. (C,D) Venn diagram and Manhattan plot for MTAG-FN-BMD. (E,F) Venn diagram and Manhattan plot for HOA; (G,H) Venn diagram and Manhattan plot for MTAG-HOA.
3.4 Genetic insights into MTAG-HOA
MAGMA analysis identified 66 pleiotropic genes surpassing the multiple testing corrected P-value threshold (Supplementary Table S3). Independently, GCTA-fastBAT pinpointed 117 pleiotropic genes meeting the same significance criterion (Supplementary Table S4). Stringent cross-validation, by taking the intersection of gene lists from both methods, revealed a high-confidence set of 48 pleiotropic genes.
3.5 GO and KEGG analysis results
Integrating GCTA-COJO, FUMA, MAGMA, and GCTA-fastBAT analyses, this study identified 13 genes as potential key genes of MTAG-HOA (Figures 4A, B; Table 2). Among these risk loci, all are situated within intronic regions of key genes, except for rs1046934, which is in an exon of TSEN15. Notably, several of these gene loci, including those within RUNX2 (rs2396502), COL11A1 (rs4338381), LTBP3 (rs10896015), SLBP (rs11732213), LMX1B (rs62578127), LRIG3 (rs79056043), MAPT (rs62063281), and COLGALT2 (rs10911472), have been previously validated in independent studies as being intricately linked to the pathogenesis of HOA.
Figure 4. Functional enrichment and tissue-specific expression analyses of key genes. (A) Identification of key genes at the intersection of genes mapped from risk loci (co-identified by GCTA-COJO and LeadSNPs) and pleiotropic genes (identified by MAGMA and GCTA-fastBAT). (B) Manhattan plot illustrating the genomic distribution of key genes associated with MTAG-HOA. (C,D) Functional enrichment analysis (GO and KEGG) of the identified key genes. (E) Tissue-specific expression of the key genes across 49 tissues, based on TWAS. (F) Tissue-specific expression of the key genes across 49 tissues, based on SMR.
Table 2. Information on key genes and risk loci.
GO analysis revealed that the key genes were significantly enriched in biological processes associated with cartilage development, osteogenesis, cell proliferation, apoptosis, and stem cell differentiation. Notably, these genes play crucial roles in the differentiation of cells that form bone, the formation of bone tissue, and the process by which bone tissue is mineralized. The analysis further highlighted regulatory mechanisms, including the TGF-β signaling pathway and pathways governing stem cell proliferation and differentiation (Figure 4C). KEGG pathway analysis demonstrated enrichment of these genes in pathways such as the MAPK signaling pathway, pathways in cancer, and inflammatory responses. This finding suggests that these genes may have pivotal roles in tumorigenesis and immune reactions (Figure 4D).
3.6 Multi -tissue TWAS identifies MTAG-HOA relevant peripheral tissues
This study employed S-PrediXcan to investigate the expression profiles of these 13 genes across 49 tissues, with the false discovery rate controlled by multiple comparison correction. Following this, Z-scores were utilized in this study to assess the strength and direction of the gene expression-phenotype association. Results revealed that ASTN2 showed no significant expression in any of the 49 tissues examined. COL11A1 downregulation occurred in five tissues, including the Adrenal Gland and Brain Putamen basal ganglia, while COLGALT2 exhibited significant upregulation in 28 tissues such as the Lung and Heart Atrial Appendage. ITIH1 downregulation was detected in the Adipose Visceral Omentum and Artery Coronary. LMX1B exhibited downregulation in the Heart Atrial Appendage and Testis, while showing upregulation in Adipose Subcutaneous and Thyroid tissues. LRIG3 displayed downregulation in seven tissues (e.g., Muscle Skeletal and Brain Anterior cingulate cortex BA24), with concurrent upregulation detected in five tissues including Breast Mammary Tissue and Minor Salivary Gland. LTBP3 demonstrated downregulation in 18 tissues (spanning the Brain Frontal Cortex BA9 and Muscle Skeletal), whereas its upregulation was specifically observed in the Testis and Pancreas. MAPT showed downregulation in the Esophagus Gastroesophageal Junction and the Small Intestine Terminal Ileum, and significant upregulation in the Spleen. RUNX2 displayed marked upregulation in the Pituitary. SLBP demonstrated reduced expression in the Minor Salivary Gland. SMAD3 downregulated in the Ovary and Muscle Skeletal, whereas upregulation was observed in the Spleen, Testis and Minor Salivary Gland. TGFA revealed downregulation in five brain regions, including the Brain Cortex. Finally, TSEN15 levels were significantly downregulated in 13 brain tissues (e.g., Brain Cortex) and 20 other non-neural tissues (Figure 4E; Supplementary Table S5).
To enhance statistical power compared to single-tissue analyses, this study employed S-MultiXcan, integrating results from multiple single-tissue models to generate unified aggregated statistics. Following multiple testing correction, 13 gene-level associations identified by S-MultiXcan were evaluated. The results showed that ASTN2 did not reach the significance threshold in the S-MultiXcan association analysis, while the remaining 12 genes demonstrated statistically significant associations (Supplementary Table S6).
3.7 SMR analysis results
SMR analysis revealed that, among the 13 genes of interest, seven genes (COLGALT2, LTBP3, MAPT, RUNX2, SMAD3, TGFA, and TSEN15) met the pre-defined selection criteria (P.adjust < 0.05 and P. HEIDI > 0.01) across 49 tissues. Specifically, COLGALT2, SMAD3, and TSEN15 demonstrated significant associations across multiple tissues. LTBP3 showed a significant association in the Pancreas tissue. MAPT showed a significant association in the Spleen tissue. RUNX2 showed a significant association in the Brain Hippocampus tissue. TGFA showed significant associations in both the Brain Caudate basal ganglia and Brain Cortex tissues (Figure 4F; Supplementary Table S7).
3.8 Proteomic studies
Plasma proteomic models from individuals of European ancestry in the INTERVAL study were pre-calculated, and PWAS analysis was conducted to identify significant associations among the 13 key genes implicated by MTAG-HOA. The analysis identified four proteins (COLGALT2, TGFA, ITIH1, and LRIG3) with significant associations. Furthermore, BLISS analysis conducted with protein data from deCODE revealed significant associations for five proteins (COLGALT2, ITIH1, LRIG3, SMAD3, and MAPT) (Table 3).
Table 3. Proteome studies.
3.9 LDSC-SEG analysis results
To investigate the potential cellular origins of MTAG-HOA, LDSC-SEG analysis was performed, referencing cell-type expression data from the Franke lab. The results indicated significant associations between MTAG-HOA and seven cell types—Cartilage, Lymphocytes, Oocytes, B Lymphocytes, Germ Cells, Osteoblasts, and Embryoid Bodies—at a coefficient level of P < 0.05 (Figure 5; Supplementary Table S8).
Figure 5. Bar plot of cell type-specific analysis using LDSC-SEG, displaying the top 50 enriched cell types ranked by p-value.
4 Discussion
4.1 Genetic correlation between HOA and FN-BMD
Emerging evidence over the past decade has implicated dynamic bone remodeling in the progression of OA, but the mechanistic links between BMD variations and joint degeneration phases remain poorly defined. Considering the pivotal role of genetic factors in both OA and its associated BMD changes, as well as their pleiotropic genetic architecture basis, traditional single disease approaches may prove insufficient to fully elucidate the underlying mechanisms. To address this gap, this study adopts a systematic and integrated strategy to comprehensively investigate the intricate genetic interactions driving the pathogenesis of both HOA and FN-BMD.
This study, employing bidirectional MR analysis, identified a negative causal effect of FN-BMD on HOA, whereas reverse MR analysis did not reveal such a causal relationship, this finding is consistent with previous MR studies (Qu et al., 2023). Furthermore, employing LDSC and HDL methods for global genetic correlation analysis, we elucidated the genetic correlation between HOA and FN-BMD, thereby providing robust evidence for their shared genetic architecture. Epidemiological studies have also indicated that osteophytic radiographic changes in HOA are associated with increased FN-BMD (Burger et al., 1996; Jadzic et al., 2022). Additionally, a longitudinal study demonstrated that radiographic changes characterized by joint space narrowing in HOA are also linked to increased BMD (Nevitt et al., 1995). The findings of this study are consistent with previous epidemiological studies, further strengthening these observations with robust genetic evidence derived from LDSC and HDL.
4.2 Risk loci, key genes and proteomics
To further increase the statistical power of the original HOA dataset, multi-trait analysis was applied to the GWAS datasets of HOA and FN-BMD. Analyses using GCTA-COJO and FUMA identified 28 independent risk loci significantly associated with MTAG-HOA, while MAGMA and GCTA-fastBAT analyses further identified 48 pleiotropic genes. Integrating the results of the GCTA-COJO, FUMA, MAGMA and GCTA-fastBAT analyses, the final refined set of 13 gene loci with significant association with MTAG-HOA was identified. These genes include COL11A1, COLGALT2, TSEN15, TGFA, ITIH1, SLBP, RUNX2, ASTN2, LMX1B, LTBP3, LRIG3, SMAD3, and MAPT. Subsequent proteomic analyses using PWAS and BLISS methods revealed key proteins significantly associated with HOA, including COLGALT2, TGFA, ITIH1, LRIG3, SMAD3, and MAPT. Integrating these findings underscores the potential importance of COLGALT2, ITIH1, and LRIG3 in HOA, and additionally reveals significant associations for TGFA, SMAD3 and MAPT. These findings not only deepen our understanding of the genetic mechanisms underlying HOA but also provide potential biomarkers and therapeutic targets for future research, thereby facilitating early diagnosis and precision treatment of HOA.
COLGALT2 encodes procollagen galactosyl transferase 2, an enzyme involved in the post-translational glycosylation of collagen (Schegg et al., 2009). Collagen is a fundamental structural and functional component of the extracellular matrix (ECM) in chondrocytes and its post-translational modifications are critical for its biological activity (Hardin et al., 2015). It has been demonstrated that increased glycosylation of collagen results in a reduction of intermolecular cross-linking, leading to a decrease in collagen fiber diameter and a concomitant reduction in tensile strength (Dominguez et al., 2005). Increased COLGALT2 expression, leading to elevated galactosyltransferase activity, is detrimental to cartilage health by impairing collagen biosynthesis (Rice et al., 2019). ITIH1 is synthesized by chondrocytes and interacts with hyaluronic acid and other extracellular matrix components, thereby conferring structural stability to cartilage (Zhao et al., 1995; Lord et al., 2020). Higher circulating ITIH1 levels have been associated with a reduced risk of HOA (Song et al., 2024). LRIG3 regulates neural crest development through modulation of the effector transcription factor Ets1 and modulates bone morphogenetic protein (BMP) signaling (Wang et al., 2015). Despite these established functions in neural development and signaling pathways, the exact role and underlying mechanisms of LRIG3 in the pathogenesis of OA remain to be thoroughly elucidated. SMAD3 maintains cartilage homeostasis and regulates inflammatory responses by activating the TGF-β signaling pathway and transducing signals from TGF-β receptors to the nucleus. Its mechanisms of action include mitochondrial dynamics, responses to mechanical stress and activation of proptosis. This signaling cascade is essential for chondrocyte function and the pathogenesis of OA (Liu et al., 2020). TGFA induces endothelin receptor A (ETA) expression in chondrocytes at both gene and protein levels, activating the RhoA/ROCK and MEK/ERK signaling pathways, ultimately driving cartilage degradation (Appleton et al., 2010). MAPT (microtubule-associated protein tau) is strongly associated with brain aging and neurodegenerative changes (Yi et al., 2025).
Although proteomic analysis revealed significant associations for the protein levels of six genes (COLGALT2, TGFA, ITIH1, LRIG3, SMAD3, and MAPT), other identified key genes did not exhibit such associations, potentially due to the limited comprehensiveness of current proteomic datasets. Consequently, despite the non-significant findings from proteomic analysis for these latter genes, their potential importance in the pathogenesis and progression of HOA should not be disregarded. RUNX2 is a transcription factor essential for skeletal system development and maturation, regulating physiological processes such as chondrocyte maturation and osteoblast differentiation (Komori, 2022). Studies have shown that RUNX2 overexpression promotes chondrocyte maturation and accelerates cartilage matrix degradation, contributing to the development of OA (Qin et al., 2015). Additionally, it facilitates the trans differentiation of terminal hypertrophic chondrocytes into osteoblasts in the femur (Takarada et al., 2013; Qin et al., 2020). Inhibition of COL11A1 results in impaired type II collagen fibril formation, subsequently affecting chondrocyte maturation and bone mineralization (Seegmiller et al., 1971). Moreover, COL11A1 suppresses lymphocyte enhancer-binding factor 1 (Lef-1), thereby inhibiting RUNX2-dependent transcriptional activation of the calcitonin promoter, ultimately affecting osteoblast differentiation and mineralization (Kahler and Westendorf, 2003). LTBP3 encodes latent TGF-β binding protein (LTBP) and plays a crucial role in TGF-β signaling regulation. Research indicates that TGF-β participates in particular chondrocyte formation and homeostatic regulation of cartilage matrix (Wang et al., 2018). Knockdown or loss of LTBP3 results in osteophyte formation, increased articular surface ossification, fibrosis, and cartilage loss, highlighting its critical role in skeletal homeostasis (Dabovic et al., 2002). LMX1B, a LIM homeobox domain transcription factor, regulates cell survival, proliferation, and inflammation. Inhibition of LMX1B enhances cell viability and suppresses chondrocyte apoptosis and inflammatory responses (Mu et al., 2022). SLBP, a cell cycle regulatory protein, facilitates the production of mature histones via 3′-end processing during the G1/S phase transition and is considered a critical gene significantly influencing joint health, particularly concerning risks associated with joint space width and hip cartilage thickness (Castaño-Betancourt et al., 2016; Wang et al., 2019). ASTN2 is highly expressed in both the developing and adult brain (Wilson et al., 2010). Studies have identified ASTN2 as a risk gene for HOA, however, its underlying mechanism remains unclear (Zeggini et al., 2012; Lindner et al., 2015). TSEN15 is a subunit of the tRNA splicing endonuclease complex. Although previous studies have indicated that TSEN15 is an independent risk gene for OA (Kehayova et al., 2023), its precise mechanisms of action remain to be further elucidated.
rs2396502, an intronic variant within the RUNX2 gene, has been identified by research as a significant risk locus for OA (Tachmazidou et al., 2019). Studies indicate that rs10896015, an intronic variant of the LTBP3 gene, is a functionally significant DNA variation significantly associated with HOA disease status at both genotype and allele levels (Zhao et al., 2020). rs10911472, an intronic variant within the COLGALT2 gene, has been suggested to influence COLGALT2 expression by modulating the methylation status of specific DNA methylation sites (cg18131582) in chondrocytes. This epigenetic regulatory mechanism may play a significant role in gene expression regulation, RNA splicing, and the generation of genetic diversity via alternative splicing (Aubourg et al., 2022). rs11732213, an intronic locus within the SLBP gene, has been demonstrated by research to have its epigenetic regulation in cartilage closely linked to SLBP function. SLBP is predominantly implicated in RNA splicing and gene expression regulation and is regarded as a pivotal gene that exerts a substantial influence on joint health, particularly about the risk associated with joint space width and hip cartilage thickness. (Castaño-Betancourt et al., 2016). The potential of rs11732213 to modulate the expression of SLBP and its associated genes by altering methylation patterns within cartilage is suggested. This regulatory activity could influence the onset and progression of HOA, thereby increasing susceptibility to the condition (Brumwell et al., 2022). rs62578127, an intronic variant within the LMX1B gene; rs4338381, an intronic variant within the COL11A1 gene; rs79056043, an intronic variant within the LRIG3 gene; and rs62063281, an intronic variant within the MAPT gene, are all considered significant risk loci for OA (Aubourg et al., 2022).
4.3 Multi-tissue analysis of potential gene expression
A multi-tissue analysis of potential gene expressions for MTAG-HOA was performed to identify and prioritize candidate causal genes and their relevant tissues. By integrating tissue-specific eQTL data from GTEx and employing TWAS (S-PrediXcan and S-MultiXcan) and SMR, the analysis identified eight gene expression regulatory loci associated with HOA across 22 tissues. Notably, brain tissue accounted for 28.6% of these associations, indicating significant enrichment. Cross-validation between transcriptome-wide association analysis and SMR highlighted potential causal associations of COLGALT2, RUNX2, and TGFA expression in brain tissue with HOA. Although the underlying mechanisms remain unclear, these findings suggest a potential role for the brain-bone axis in HOA pathogenesis, aligning with recent advancements in this emerging field. The brain-bone axis represents a complex bidirectional communication network linking the central nervous system and bone metabolism through diverse signaling pathways, neuroendocrine factors, and molecular mediators (Quiros-Gonzalez and Yadav, 2014). Growing evidence indicates that multiple neural pathways, including the sympathetic nervous system, hypothalamic neuropeptides, and neurotransmitters, play a crucial role in regulating both physiological and pathological bone changes (Karsenty and Ferron, 2012).
4.4 Enrichment analysis and cell specificity
This study employed GO and KEGG enrichment analyses, alongside LDSC-SEG cell-type specificity analysis, to thoroughly investigate the genetic relationship between HOA and FN-BMD and their underlying biological mechanisms. GO and KEGG enrichment analyses revealed that key genes were significantly enriched in pathways associated with biological processes crucial for skeletal development, cartilage formation, cell proliferation, apoptosis, and stem cell differentiation. These findings underscore the critical involvement of these genes in pivotal processes such as osteogenesis, bone tissue mineralization, and osteocyte differentiation, particularly within the TGF-β signaling pathway and mechanisms governing stem cell proliferation and differentiation. Concurrently, LDSC-SEG cell-type specificity analysis provided further insights into the cellular origins of HOA-associated genes. Integration with cell type-specific expression data revealed significant associations of MTAG-identified HOA loci with multiple cell types, including chondrocytes, osteoblasts, and immune cells such as lymphocytes and B lymphocytes. The associations with chondrocytes and osteoblasts particularly underscore the central role of these genes in bone tissue formation, whereas the relevance to immune cells suggests their potential pivotal involvement in immune responses and the inflammatory processes of osteoarthritis. These discoveries offer novel perspectives on the pathogenesis of HOA, particularly concerning how these genes function within different cell types and their interplay between the skeletal and immune systems.
4.5 Advantages and limitations
In summary, this study provides genetic evidence linking HOA and FN-BMD, while unraveling the genetic landscape of HOA. The strengths of this study are manifold. First, it integrates LDSC and HDL, leveraging the unique strengths of each method to minimize false-negative and false-positive findings. Second, the use of MTAG for multitrait GWAS analysis on the HOA and FN-BMD pathway represents a pioneering effort in this domain. Third, GCTA-COJO and FUMA revealed independent risk loci for MTAG-HOA, while MAGMA and GCTA-fastBAT identified 13 pleiotropic genes. Fourth, TWAS and SMR analyses identified eight genes whose expression regulation was associated with HOA risk across 22 tissues. Of note, brain tissue accounted for a substantial proportion of these associations.
It is imperative to acknowledge the limitations of this study. First, the focus on European populations necessitates cautious interpretation when extending conclusions to diverse racial and ethnic groups. Second, although TWAS and SMR analyses were conducted using the GTEx V8 eQTL dataset, only 49 tissues were included, and key tissues such as synovial tissue were not covered. Finally, although MTAG-HOA pleiotropic genes were identified, further longitudinal and experimental studies are crucial for elucidating the underlying biological mechanisms.
5 Conclusion
This study elucidates the genetic link between HOA and FN-BMD, highlights pleiotropic genes in HOA, and uncovers novel genetic risk factors. These findings deepen our understanding of HOA’s genetic landscape and may pave the way for developing innovative therapeutic strategies.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Author contributions
JZ: Conceptualization, Data curation, Formal Analysis, Project administration, Writing – review and editing. JN: Formal analysis, Visualization, Data curation, Writing – review and editing. ZJ: Conceptualization, Data curation, Formal Analysis, Writing – review and editing. XD: Conceptualization, Data curation, Project administration, Resources, Writing – review and editing. SW: Conceptualization, Data curation, Formal Analysis, Project administration, Investigation, Methodology, Supervision, Writing – review and editing. HL: Conceptualization, Data curation, Formal Analysis, Methodology, Project administration, Supervision, Funding acquisition, Writing – review and editing. JK: Conceptualization, Funding acquisition, Supervision, Writing -review and editing.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. National Famous Liaoning Provincial Science and Technology Programme Joint Programme (Natural Science Foundation - top level project) (2024-MSLH-261).
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 author(s) 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/fgene.2025.1597005/full#supplementary-material
References
Appleton, C. T., Usmani, S. E., Mort, J. S., and Beier, F. (2010). Rho/ROCK and MEK/ERK activation by transforming growth factor-alpha induces articular cartilage degradation. Lab. Invest 90 (1), 20–30. doi:10.1038/labinvest.2009.111
Aubourg, G., Rice, S. J., Bruce-Wootton, P., and Loughlin, J. (2022). Genetics of osteoarthritis. Osteoarthr. Cartil. 30 (5), 636–649. doi:10.1016/j.joca.2021.03.002
Auger, J. D., Naik, A. J., Murakami, A. M., Gerstenfeld, L. C., and Morgan, E. F. (2022). Spatial assessment of femoral neck bone density and microstructure in hip osteoarthritis. Bone Rep. 16, 101155. doi:10.1016/j.bonr.2021.101155
Auton, A., Brooks, L. D., Durbin, R. M., Garrison, E. P., Kang, H. M., Korbel, J. O., et al. (2015). A global reference for human genetic variation. Nature 526 (7571), 68–74. doi:10.1038/nature15393
Bakshi, A., Zhu, Z., Vinkhuyzen, A. A., Hill, W. D., McRae, A. F., Visscher, P. M., et al. (2016). Fast set-based association analysis using summary data from GWAS identifies novel gene loci for human complex traits. Sci. Rep. 6, 32894. doi:10.1038/srep32894
Barbeira, A. N., Dickinson, S. P., Bonazzola, R., Zheng, J., Wheeler, H. E., Torres, J. M., et al. (2018). Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics. Nat. Commun. 9 (1), 1825. doi:10.1038/s41467-018-03621-1
Barbeira, A. N., Pividori, M., Zheng, J., Wheeler, H. E., Nicolae, D. L., and Im, H. K. (2019). Integrating predicted transcriptome from multiple tissues improves association detection. PLoS Genet. 15 (1), e1007889. doi:10.1371/journal.pgen.1007889
Brumwell, A., Aubourg, G., Hussain, J., Parker, E., Deehan, D. J., Rice, S. J., et al. (2022). Identification of TMEM129, encoding a ubiquitin-protein ligase, as an effector gene of osteoarthritis genetic risk. Arthritis Res. Ther. 24 (1), 189. doi:10.1186/s13075-022-02882-y
Bulik-Sullivan, B., Finucane, H. K., Anttila, V., Gusev, A., Day, F. R., Loh, P. R., et al. (2015). An atlas of genetic correlations across human diseases and traits. Nat. Genet. 47 (11), 1236–1241. doi:10.1038/ng.3406
Burger, H., van Daele, P. L., Odding, E., Valkenburg, H. A., Hofman, A., Grobbee, D. E., et al. (1996). Association of radiographically evident osteoarthritis with higher bone mineral density and increased bone loss with age. The rotterdam study. Arthritis Rheum. 39 (1), 81–86. doi:10.1002/art.1780390111
Carpenter, M. R., Carpenter, R. L., Peel, J., Zukley, L. M., Angelopoulou, K. M., Fischer, I., et al. (2006). The reliability of isokinetic and isometric leg strength measures among individuals with symptoms of mild osteoarthritis. J. Sports Med. Phys. Fit. 46 (4), 585–589.
Castaño-Betancourt, M. C., Evans, D. S., Ramos, Y. F., Boer, C. G., Metrustry, S., Liu, Y., et al. (2016). Novel genetic variants for cartilage thickness and hip osteoarthritis. PLoS Genet. 12 (10), e1006260. doi:10.1371/journal.pgen.1006260
Chiba, K., Burghardt, A. J., Osaki, M., and Majumdar, S. (2014). Three-dimensional analysis of subchondral cysts in hip osteoarthritis: an ex vivo HR-pQCT study. Bone 66, 140–145. doi:10.1016/j.bone.2014.06.001
Cui, A., Li, H., Wang, D., Zhong, J., Chen, Y., and Lu, H. (2020). Global, regional prevalence, incidence and risk factors of knee osteoarthritis in population-based studies. EClinicalMedicine 29-30, 100587. doi:10.1016/j.eclinm.2020.100587
Dabovic, B., Chen, Y., Colarossi, C., Obata, H., Zambuto, L., Perle, M. A., et al. (2002). Bone abnormalities in latent TGF-[beta] binding protein (Ltbp)-3-null mice indicate a role for Ltbp-3 in modulating TGF-[beta] bioavailability. J. Cell Biol. 156 (2), 227–232. doi:10.1083/jcb.200111080
de Leeuw, C. A., Mooij, J. M., Heskes, T., and Posthuma, D. (2015). MAGMA: generalized gene-set analysis of GWAS data. PLoS Comput. Biol. 11 (4), e1004219. doi:10.1371/journal.pcbi.1004219
Dominguez, L. J., Barbagallo, M., and Moro, L. (2005). Collagen overglycosylation: a biochemical feature that May contribute to bone quality. Biochem. Biophys. Res. Commun. 330 (1), 1–4. doi:10.1016/j.bbrc.2005.02.050
Finucane, H. K., Reshef, Y. A., Anttila, V., Slowikowski, K., Gusev, A., Byrnes, A., et al. (2018). Heritability enrichment of specifically expressed genes identifies disease-relevant tissues and cell types. Nat. Genet. 50 (4), 621–629. doi:10.1038/s41588-018-0081-4
Guo, P., Gong, W., Li, Y., Liu, L., Yan, R., Wang, Y., et al. (2022). Pinpointing novel risk loci for lewy body dementia and the shared genetic etiology with alzheimer's disease and parkinson's disease: a large-scale multi-trait association analysis. BMC Med. 20 (1), 214. doi:10.1186/s12916-022-02404-2
Hardin, J. A., Cobelli, N., and Santambrogio, L. (2015). Consequences of metabolic and oxidative modifications of cartilage tissue. Nat. Rev. Rheumatol. 11 (9), 521–529. doi:10.1038/nrrheum.2015.70
Husted, H. (2012). Fast-track hip and knee arthroplasty: clinical and organizational aspects. Acta Orthop. Suppl. 83 (346), 1–39. doi:10.3109/17453674.2012.700593
Jadzic, J., Zagorac, S., Djuric, M., and Djonic, D. (2022). Hip structure analysis and femoral osteodensitometry in aged postmenopausal women with hip osteoarthritis and femoral neck fracture. Int. Orthop. 46 (12), 2747–2755. doi:10.1007/s00264-022-05602-y
Kahler, R. A., and Westendorf, J. J. (2003). Lymphoid enhancer factor-1 and beta-catenin inhibit Runx2-dependent transcriptional activation of the osteocalcin promoter. J. Biol. Chem. 278 (14), 11937–11944. doi:10.1074/jbc.M211443200
Karsenty, G., and Ferron, M. (2012). The contribution of bone to whole-organism physiology. Nature 481 (7381), 314–320. doi:10.1038/nature10763
Kehayova, Y. S., Wilkinson, J. M., Rice, S. J., and Loughlin, J. (2023). Mediation of the same epigenetic and transcriptional effect by independent osteoarthritis risk-conferring alleles on a shared target gene, COLGALT2. Arthritis Rheumatol. 75 (6), 910–922. doi:10.1002/art.42427
Komori, T. (2022). Whole aspect of Runx2 functions in skeletal development. Int. J. Mol. Sci. 23 (10), 5776. doi:10.3390/ijms23105776
Lee, C. H., Shi, H., Pasaniuc, B., Eskin, E., and Han, B. (2021). PLEIO: a method to map and interpret pleiotropic loci with GWAS summary statistics. Am. J. Hum. Genet. 108 (1), 36–48. doi:10.1016/j.ajhg.2020.11.017
Li, Y. I., Knowles, D. A., Humphrey, J., Barbeira, A. N., Dickinson, S. P., Im, H. K., et al. (2018). Annotation-free quantification of RNA splicing using LeafCutter. Nat. Genet. 50 (1), 151–158. doi:10.1038/s41588-017-0004-9
Lin, S. H., Brown, D. W., and Machiela, M. J. (2020). LDtrait: an online tool for identifying published phenotype associations in linkage disequilibrium. Cancer Res. 80 (16), 3443–3446. doi:10.1158/0008-5472.Can-20-0985
Lindner, C., Thiagarajah, S., Wilkinson, J. M., Panoutsopoulou, K., Day-Williams, A. G., Cootes, T. F., et al. (2015). Investigation of association between hip osteoarthritis susceptibility loci and radiographic proximal femur shape. Arthritis Rheumatol. 67 (8), 2076–2084. doi:10.1002/art.39186
Liu, J., Chen, Q., Alkam, E., Zheng, X., Li, Y., Wang, L., et al. (2020). Association between gene polymorphisms of TGF-β and Smad3 and susceptibility to arthritis: a meta-analysis. Expert Rev. Clin. Immunol. 16 (9), 943–954. doi:10.1080/1744666x.2020.1816826
Liu, L., Tian, F., Li, G. Y., Xu, W., and Xia, R. (2022). The effects and significance of gut microbiota and its metabolites on the regulation of osteoarthritis: close coordination of gut-bone axis. Front. Nutr. 9, 1012087. doi:10.3389/fnut.2022.1012087
Liu, X., Li, D., Gao, W., Liu, H., Chen, P., Zhao, Y., et al. (2024). Shared genetic architecture between COVID-19 and irritable bowel syndrome: a large-scale genome-wide cross-trait analysis. Front. Immunol. 15, 1442693. doi:10.3389/fimmu.2024.1442693
Lord, M. S., Melrose, J., Day, A. J., and Whitelock, J. M. (2020). The Inter-α-Trypsin inhibitor family: versatile molecules in biology and pathology. J. Histochem Cytochem 68 (12), 907–927. doi:10.1369/0022155420940067
Meng, Z., Zhang, C., Liu, S., Li, W., Wang, Y., Zhang, Q., et al. (2025). Exploring genetic loci linked to COVID-19 severity and immune response through multi-trait GWAS analyses. Front. Genet. 16, 1502839. doi:10.3389/fgene.2025.1502839
Motta, F., Barone, E., Sica, A., and Selmi, C. (2023). Inflammaging and osteoarthritis. Clin. Rev. Allergy Immunol. 64 (2), 222–238. doi:10.1007/s12016-022-08941-1
Mu, Y., Wang, L., Fu, L., and Li, Q. (2022). Knockdown of LMX1B suppressed cell apoptosis and inflammatory response in IL-1β-Induced human osteoarthritis chondrocytes through NF-κB and NLRP3 signal pathway. Mediat. Inflamm. 2022, 1870579. doi:10.1155/2022/1870579
Nevitt, M. C., Lane, N. E., Scott, J. C., Hochberg, M. C., Pressman, A. R., Genant, H. K., et al. (1995). Radiographic osteoarthritis of the hip and bone mineral density. The study of osteoporotic fractures research group. Arthritis Rheum. 38 (7), 907–916. doi:10.1002/art.1780380706
Ning, Z., Pawitan, Y., and Shen, X. (2020). High-definition likelihood inference of genetic correlations across human complex traits. Nat. Genet. 52 (8), 859–864. doi:10.1038/s41588-020-0653-y
Palmer, T. M., Lawlor, D. A., Harbord, R. M., Sheehan, N. A., Tobias, J. H., Timpson, N. J., et al. (2012). Using multiple genetic variants as instrumental variables for modifiable risk factors. Stat. Methods Med. Res. 21 (3), 223–242. doi:10.1177/0962280210394459
Qin, X., Jiang, Q., Matsuo, Y., Kawane, T., Komori, H., Moriishi, T., et al. (2015). Cbfb regulates bone development by stabilizing runx family proteins. J. Bone Min. Res. 30 (4), 706–714. doi:10.1002/jbmr.2379
Qin, X., Jiang, Q., Nagano, K., Moriishi, T., Miyazaki, T., Komori, H., et al. (2020). Runx2 is essential for the transdifferentiation of chondrocytes into osteoblasts. PLoS Genet. 16 (11), e1009169. doi:10.1371/journal.pgen.1009169
Qu, Y., Chen, S., Han, M., Gu, Z., Zhang, Y., Fan, T., et al. (2023). Osteoporosis and osteoarthritis: a bi-directional Mendelian randomization study. Arthritis Res. Ther. 25 (1), 242. doi:10.1186/s13075-023-03213-5
Quiros-Gonzalez, I., and Yadav, V. K. (2014). Central genes, pathways and modules that regulate bone mass. Arch. Biochem. Biophys. 561, 130–136. doi:10.1016/j.abb.2014.06.005
Rice, S. J., Cheung, K., Reynard, L. N., and Loughlin, J. (2019). Discovery and analysis of methylation quantitative trait loci (mQTLs) mapping to novel osteoarthritis genetic risk signals. Osteoarthr. Cartil. 27 (10), 1545–1556. doi:10.1016/j.joca.2019.05.017
Schegg, B., Hülsmeier, A. J., Rutschmann, C., Maag, C., and Hennet, T. (2009). Core glycosylation of collagen is initiated by two beta(1-O)galactosyltransferases. Mol. Cell Biol. 29 (4), 943–952. doi:10.1128/mcb.02085-07
Seegmiller, R., Fraser, F. C., and Sheldon, H. (1971). A new chondrodystrophic mutant in mice. Electron microscopy of normal and abnormal chondrogenesis. J. Cell Biol. 48 (3), 580–593. doi:10.1083/jcb.48.3.580
Shan, Y., Hu, H., and Chu, Y. (2025). Cross-ancestry genome-wide association study identifies new susceptibility genes for preeclampsia. BMC Pregnancy Childbirth 25 (1), 379. doi:10.1186/s12884-025-07534-y
Sivakumaran, S., Agakov, F., Theodoratou, E., Prendergast, J. G., Zgaga, L., Manolio, T., et al. (2011). Abundant pleiotropy in human complex diseases and traits. Am. J. Hum. Genet. 89 (5), 607–618. doi:10.1016/j.ajhg.2011.10.004
Song, S., Qiao, J., Zhao, R., Lu, Y. J., Wang, C., Chang, M. J., et al. (2024). Identification of novel drug targets for osteoarthritis by integrating genetics and proteomes from blood. J. Orthop. Surg. Res. 19 (1), 559. doi:10.1186/s13018-024-05034-x
Song, Y., Li, L., Jiang, Y., Peng, B., Jiang, H., Chao, Z., et al. (2025). Multitrait genetic analysis identifies novel pleiotropic loci for depression and schizophrenia in east Asians. Schizophr. Bull. 51 (3), 684–695. doi:10.1093/schbul/sbae145
Sun, B. B., Maranville, J. C., Peters, J. E., Stacey, D., Staley, J. R., Blackshaw, J., et al. (2018). Genomic atlas of the human plasma proteome. Nature 558 (7708), 73–79. doi:10.1038/s41586-018-0175-2
Tachmazidou, I., Hatzikotoulas, K., Southam, L., Esparza-Gordillo, J., Haberland, V., Zheng, J., et al. (2019). Identification of new therapeutic targets for osteoarthritis through genome-wide analyses of UK biobank data. Nat. Genet. 51 (2), 230–236. doi:10.1038/s41588-018-0327-1
Takarada, T., Hinoi, E., Nakazato, R., Ochi, H., Xu, C., Tsuchikane, A., et al. (2013). An analysis of skeletal development in osteoblast-specific and chondrocyte-specific runt-related transcription factor-2 (Runx2) knockout mice. J. Bone Min. Res. 28 (10), 2064–2069. doi:10.1002/jbmr.1945
Turley, P., Walters, R. K., Maghzian, O., Okbay, A., Lee, J. J., Fontana, M. A., et al. (2018). Multi-trait analysis of genome-wide association summary statistics using MTAG. Nat. Genet. 50 (2), 229–237. doi:10.1038/s41588-017-0009-4
Urbut, S. M., Wang, G., Carbonetto, P., and Stephens, M. (2019). Flexible statistical methods for estimating and testing effects in genomic studies with multiple conditions. Nat. Genet. 51 (1), 187–195. doi:10.1038/s41588-018-0268-8
Verbanck, M., Chen, C. Y., Neale, B., and Do, R. (2018). Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nat. Genet. 50 (5), 693–698. doi:10.1038/s41588-018-0099-7
Wang, C., Kam, R. K., Shi, W., Xia, Y., Chen, X., Cao, Y., et al. (2015). The proto-oncogene transcription factor Ets1 regulates neural crest development through histone deacetylase 1 to mediate output of bone morphogenetic protein signaling. J. Biol. Chem. 290 (36), 21925–21938. doi:10.1074/jbc.M115.644864
Wang, B., Overgaard, S., Chemnitz, J., and Ding, M. (2016). Cancellous and cortical bone microarchitectures of femoral neck in rheumatoid arthritis and osteoarthritis compared with donor controls. Calcif. Tissue Int. 98 (5), 456–464. doi:10.1007/s00223-015-0098-y
Wang, C., Silverman, R. M., Shen, J., and O'Keefe, R. J. (2018). Distinct metabolic programs induced by TGF-β1 and BMP2 in human articular chondrocytes with osteoarthritis. J. Orthop. Transl. 12, 66–73. doi:10.1016/j.jot.2017.12.004
Wang, Y., Huang, H., and Li, Y. (2019). Knocking Down miR-384 promotes growth and metastasis of osteosarcoma MG63 cells by targeting SLBP. Artif. Cells Nanomed Biotechnol. 47 (1), 1458–1465. doi:10.1080/21691401.2019.1601099
Watanabe, K., Taskesen, E., van Bochoven, A., and Posthuma, D. (2017). Functional mapping and annotation of genetic associations with FUMA. Nat. Commun. 8 (1), 1826. doi:10.1038/s41467-017-01261-5
Wen, X., Lee, Y., Luca, F., and Pique-Regi, R. (2016). Efficient integrative Multi-SNP association analysis via deterministic approximation of posteriors. Am. J. Hum. Genet. 98 (6), 1114–1129. doi:10.1016/j.ajhg.2016.03.029
Wilson, P. M., Fryer, R. H., Fang, Y., and Hatten, M. E. (2010). Astn2, a novel member of the astrotactin gene family, regulates the trafficking of ASTN1 during glial-guided neuronal migration. J. Neurosci. 30 (25), 8529–8540. doi:10.1523/jneurosci.0032-10.2010
Xu, C., Wang, S., Chen, X., Zhang, T., Ni, Z., Ji, G., et al. (2024). Causal associations between circulating immune cells and osteoarthritis: a bidirectional mendelian randomization study. Int. Immunopharmacol. 142 (Pt A), 113156. doi:10.1016/j.intimp.2024.113156
Yang, J., Ferreira, T., Morris, A. P., Medland, S. E., Madden, P. A., Heath, A. C., et al. (2012). Conditional and joint multiple-SNP analysis of GWAS summary statistics identifies additional variants influencing complex traits. Nat. Genet. 44 (4), 369–S3. doi:10.1038/ng.2213
Yi, F., Yuan, J., Somekh, J., Peleg, M., Zhu, Y. C., Jia, Z., et al. (2025). Genetically supported targets and drug repurposing for brain aging: a systematic study in the UK biobank. Sci. Adv. 11 (11), eadr3757. doi:10.1126/sciadv.adr3757
Yoshino, K., Hagiwara, S., Nakamura, J., Horikoshi, T., Yokota, H., Shimokawa, K., et al. (2021). The distance between the femoral nerve and anterior acetabulum is significantly shorter in hip osteoarthritis than in non-osteoarthritis hip. BMC Musculoskelet. Disord. 22 (1), 416. doi:10.1186/s12891-021-04295-5
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16 (5), 284–287. doi:10.1089/omi.2011.0118
Zapata-Linares, N., Toillon, I., Wanherdrick, K., Pigenet, A., Duhalde, F., Binvignat, M., et al. (2025). Implication of bone marrow adipose tissue in bone homeostasis during osteoarthritis. Osteoarthr. Cartil. 33, 951–964. doi:10.1016/j.joca.2025.03.004
Zeggini, E., Panoutsopoulou, K., Southam, L., Rayner, N. W., Day-Williams, A. G., Lopes, M. C., et al. (2012). Identification of new susceptibility loci for osteoarthritis (arcOGEN): a genome-wide association study. Lancet 380 (9844), 815–823. doi:10.1016/s0140-6736(12)60681-3
Zengini, E., Hatzikotoulas, K., Tachmazidou, I., Steinberg, J., Hartwig, F. P., Southam, L., et al. (2018). Genome-wide analyses using UK biobank data provide insights into the genetic architecture of osteoarthritis. Nat. Genet. 50 (4), 549–558. doi:10.1038/s41588-018-0079-y
Zhang, J., Li, Y., Li, Y., and Liu, H. (2025). Unraveling the brain-joint axis: genetic, transcriptomic, and cohort insights from neuroticism to osteoarthritis. Mamm. Genome 36, 638–650. doi:10.1007/s00335-025-10112-4
Zhao, M., Yoneda, M., Ohashi, Y., Kurono, S., Iwata, H., Ohnuki, Y., et al. (1995). Evidence for the covalent binding of SHAP, heavy chains of inter-alpha-trypsin inhibitor, to hyaluronan. J. Biol. Chem. 270 (44), 26657–26663. doi:10.1074/jbc.270.44.26657
Zhao, T., Zhao, J., Ma, C., Wei, J., Wei, B., and Liu, J. (2020). Common variants in LTBP3 gene contributed to the risk of hip osteoarthritis in han Chinese population. Biosci. Rep. 40 (6). doi:10.1042/bsr20192999
Zheng, H. F., Forgetta, V., Hsu, Y. H., Estrada, K., Rosello-Diez, A., Leo, P. J., et al. (2015). Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture. Nature 526 (7571), 112–117. doi:10.1038/nature14878
Zhu, Z., Zhang, F., Hu, H., Bakshi, A., Robinson, M. R., Powell, J. E., et al. (2016). Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets. Nat. Genet. 48 (5), 481–487. doi:10.1038/ng.3538
Glossary
OA Osteoarthritis
HOA Hip Osteoarthritis
FN-BMD Femoral Neck Bone Mineral Density
MR Mendelian randomization
LDSC Linkage Disequilibrium Score Regression
HDL High-Definition Likelihood Inference
GWAS Genome-Wide Association Studie
GCTA-COJO conditional and joint analysis
FUMA Functional Mapping and Annotation
GO Gene Ontology
KEGG Kyoto Encyclopedia of Genes and Genomes
TWAS Transcriptome-Wide Association study
BMD Bone Mineral Density
GEFOS Genetic Factors for Osteoporosis Consortium
PWAS Proteome-Wide Association Study
BLISS Biomarker Level Imputation from Summary Statistics
SNP Single-Nucleotide Polymorphism
se Standard Error
MTAG Multi-trait analysis of GWAS
DXA dual-energy X-ray absorptiometry
MAF minor allele frequency
BMI body mass index
IV instrumental variables
IVW Inverse Variance Weighted
BP biological processes
CC cellular components
MF molecular functions
BH Benjamini-Hochberg
FDR False Discovery Rate
SMR summary-data-based MR
S-LDSC Stratified LD Score Regression
Lef-1 Lymphocyte Enhancer-Binding Factor 1
ECM Extracellular Matrix
ETA Endothelin Receptor A
BMP Bone Morphogenetic Protein
Keywords: hip osteoarthritis, femoral neck bone mineral density, genetic loci, genetic heterogeneity, multi-tissue transcriptome-wide association studies
Citation: Zhou J, Na J, Jiang Z, Dou X, Wang S, Li H and Kang J (2025) Deciphering the shared genetic structure between hip osteoarthritis and femoral neck bone mineral density. Front. Genet. 16:1597005. doi: 10.3389/fgene.2025.1597005
Received: 20 March 2025; Accepted: 28 August 2025;
Published: 19 September 2025.
Edited by:
Jared C. Roach, Institute for Systems Biology (ISB), United StatesReviewed by:
Jiawen Xu, Johns Hopkins University, United StatesEucharist Kun, The University of Texas at Austin, United States
Copyright © 2025 Zhou, Na, Jiang, Dou, Wang, Li and Kang. 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: Shixuan Wang, d2FuZ3NoaXh1YW45OUAxMjYuY29t; Hongtao Li, OTQyMTMzNzJAcXEuY29t; Jian Kang, a2ozMTY1MDMwNDRAZ21haWwuY29t
†These authors have contributed equally to this work and share first authorship