Integrative Transcriptome-Wide Association Study and Mendelian Randomization Identify LYNX1 and MS4A14 as Therapeutic Targets for Osteomyelitis

Summary

This study integrates transcriptome-wide association analysis and Mendelian randomization to systematically identify drug targets for osteomyelitis. Leveraging genetic databases, the authors identified LYNX1 and MS4A14 as key candidate genes, highlighting critical molecular pathways and informing precise therapeutic strategies for managing bone infections.

Abstract

Osteomyelitis is a severe infectious disease characterized by profound inflammation of the bone and bone marrow and is predominantly caused by bacterial pathogens such as Staphylococcus aureus (S. aureus). However, current treatment regimens are significantly hindered by the emergence of antibiotic resistance and recurrent infections, driven by robust biofilm formation and intracellular bacterial persistence. Thus, there is an urgent need to identify novel therapeutic targets beyond conventional antimicrobials. An integrative computational approach combining a transcriptome-wide association study (TWAS) with summary-based Mendelian randomization (SMR) analysis was applied to systematically identify causal susceptibility genes underlying osteomyelitis. Large-scale genetic variants were leveraged as instrumental variables to predict gene expression across tissues, enabling the exploration of causal associations between these targets and disease risk. The integrative framework identified LYNX1 and MS4A14 as key candidate genes and potential therapeutic targets. Specifically, LYNX1 was associated with increased susceptibility to osteomyelitis, whereas MS4A14 exhibited potential protective properties. These findings highlight the regulatory roles of these genes in host immune response and inflammatory modulation during bone infection. This study bridges the gap between genome-wide association findings and biological interpretation, advancing the understanding of the genetic basis of osteomyelitis and supporting the development of targeted precision medicine strategies to enhance host defense and overcome therapeutic resistance.

Introduction

Osteomyelitis is a severe infectious disease characterized by profound inflammation of the bone and bone marrow. It is frequently initiated by bacterial invasion, with Staphylococcus aureus (S. aureus) being a common pathogen across acute and chronic hematogenous presentations^1,2. While S. aureus is prominent, osteomyelitis is often polymicrobial, and eradication of S. aureus alone does not always resolve the infection due to the persistence of uncultivable microorganisms and dynamic biofilm interactions^3,4. In susceptible hosts, pathogens employ intricate evasion strategies, including intracellular persistence within macrophages and the formation of dense biofilms. These biofilms, composed of extracellular DNA, polysaccharides, and specialized proteins, create formidable physical barriers that resist host immune defenses and limit the penetration of conventional antimicrobial agents^5,6. Consequently, although clinical management typically involves aggressive surgical debridement and prolonged antibiotic therapy, the presence of resilient biofilms and polymicrobial communities significantly reduces treatment efficacy, leading to prolonged morbidity and recurrent infections^7,8. Elucidation of the molecular architecture underlying host–pathogen interactions remains critical for identifying novel therapeutic targets⁹.

Genome-wide association studies (GWAS) have advanced the understanding of osteomyelitis by identifying multiple genetic loci associated with disease susceptibility¹⁰. However, due to stringent multiple testing thresholds and limitations in statistical power and sample size, many critical loci influencing osteomyelitis susceptibility are likely to remain undetected⁷. Furthermore, interpretation of GWAS findings remains challenging, as most identified variants are located in non-coding or intergenic regions, suggesting regulatory effects on gene expression rather than direct alterations in protein structure¹¹. In addition, primary GWAS data lack direct tissue specificity, necessitating complementary downstream analyses to determine the cellular contexts in which these variants exert their effects. Transcriptome-wide association studies (TWAS) address these limitations by integrating genetically predicted gene expression from reference panels into existing GWAS datasets, enabling systematic identification of functionally relevant disease-associated genes across multiple tissues¹².

To strengthen causal inference, Mendelian randomization (MR) employs genetic variants as instrumental variables, leveraging the random allocation of alleles at conception to reduce confounding and reverse causation^13,14. Expression quantitative trait loci (eQTL) analyses further link GWAS variants to gene transcription¹⁵, while summary-based Mendelian randomization (SMR) integrates GWAS and QTL datasets to prioritize causal genes and distinguish them from associations driven by linkage disequilibrium through heterogeneity testing¹⁶.

This study represents the first application of an integrative TWAS and MR framework specifically focused on osteomyelitis, addressing the gap between GWAS findings and clinically actionable insights. By combining experimental transcriptomic data with large-scale population-genetic data, disease-susceptibility genes were systematically mapped and causal relationships inferred. This integrative strategy aims to uncover previously unrecognized molecular pathways involved in host resistance to osteomyelitis and to identify candidate targets for precision therapeutic development.

Protocol

All summary statistics utilized in the Mendelian Randomization (MR) and Transcriptome-Wide Association Study (TWAS) analyses were derived strictly from previously published, de-identified datasets. Ethical approval and individual consent for the original studies are documented in their respective publications. Consequently, additional ethical approval for this data-mining study was waived by the Institutional Review Board of Tongde Hospital of Zhejiang Province (Zhe Tongde Lunshen 2024 [Yan] No. 028-JY). The tools used for this research are listed in the Table of Materials.

1. RNA-seq data acquisition and processing

Transcriptome data were obtained from the Gene Expression Omnibus (GEO) database (GSE272198) to assess conservation of innate immune pathways across mammalian species for initial validation¹⁷. Bone marrow-derived macrophages (BMDMs) were infected with S. aureus (multiplicity of infection, MOI = 10) for 1 h, followed by treatment with lysostaphin (20 µg/mL) and gentamicin (50 µg/mL) to remove extracellular bacteria. After three washes with phosphate-buffered saline (PBS), BMDMs were cultured for 24 h, lysed in a total RNA extraction reagent, and sequenced.

RNA quality was assessed using an automated electrophoresis system to ensure integrity. Libraries were prepared from three independent experiments and sequenced on a high-throughput sequencing platform. Raw reads were aligned to the mouse genome (GRCm38, mm10) using STAR (v2.7.10a). Differentially expressed genes (DEGs) were identified using DESeq2 (v1.38.0). To mitigate false positives, statistical significance was defined as an adjusted p-value (FDR) < 0.05 and |log₂ fold change| > 1. Gene Ontology (GO) analysis was performed using clusterProfiler (v4.6.0), and Gene Set Enrichment Analysis (GSEA) was conducted using GseaVis (v0.0.5). Heatmaps were generated using the pheatmap package (v1.0.12) in R (v4.2.0).

TWAS analysis

Whole-blood RNA sequencing and whole-genome sequencing (WGS) data were obtained from the Genotype-Tissue Expression (GTEx) project (V8)¹⁸. Pre-trained gene expression models were utilized from a public repository (https://doi.org/10.5281/zenodo.3842289). Osteomyelitis summary statistics for TWAS were retrieved from the FinnGen consortium, comprising 2,336 cases and 473,264 controls¹².

TWAS was conducted using three algorithms: joint-tissue imputation (JTI), PrediXcan¹⁹, and UTMOST^12,20. JTI estimates gene expression similarity and epigenetic chromatin accessibility to optimize prediction accuracy. PrediXcan applies elastic net regression with fivefold cross-validation, while UTMOST enhances accuracy by leveraging multi-tissue expression data using sparse group LASSO. The modified UTMOST framework described by Zhou et al.¹² standardizes hyperparameters for unbiased estimation. Genes with stable cross-validation scores—pre-defined as a correlation coefficient r > 0.1 and predictive significance p < 0.05²¹—were retained as imputable. Whole-blood transcriptome models were established using SNP covariance matrices from the 1000 Genomes reference dataset.

Associations between predicted gene expression and osteomyelitis risk were subsequently analyzed. To account for multiple testing, statistical significance for TWAS was primarily defined using a False Discovery Rate (FDR) threshold of < 0.05. Given the hypothesis-generating nature of this multi-stage study, loci meeting a suggestive (nominal) threshold of p < 0.05 were also prioritized for downstream Mendelian randomization (SMR) and colocalization analyses. This integrative strategy aims to maximize the capture of potential regulatory drivers while relying on multi-omic cross-validation (TWAS + SMR) to ensure the robustness of the prioritized candidates.

SMR Analysis

This study adhered to the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) guidelines²². To computationally define a phenotype representing genetic predisposition to mitochondrial dysfunction (hereafter termed “mitodys” for analysis purposes), transcripts corresponding to all known mitochondrial-related genes were extracted from the MitoCarta3.0 database²³. This gene set served as a predefined, biology-informed basis for subsequent polygenic risk prediction. All downstream functional interpretations relating to “mitodys” are derived from this computational inference and should be considered predictive and hypothesis-generating.

Expression quantitative trait loci (eQTL) instruments were generated using variants within 1000 kb of coding sequences (cis-eQTLs). Summary statistics were sourced from the eQTLGen Consortium and GTEx V8²⁴. A total of 8,932,843 SNPs linked to 1,013 mitodys-related transcripts were selected based on a genome-wide significance threshold P < 5E-8. Baseline GWAS statistics for osteomyelitis outcomes were obtained from FinnGen²⁰.

Summary-data-based Mendelian Randomization (SMR) analysis was performed using SMR (version 1.0.3) with default parameters to estimate pleiotropic associations between gene expression traits and osteomyelitis outcomes. The causal effect beta_mitodys–osteomyelitis represents the estimated log-odds effect size of mitochondrial dysfunction on osteomyelitis and is calculated as:

Odds ratios (ORs) represent the change per one-unit natural logarithmic increase in standardized gene expression levels. Co-localization was further evaluated using the heterogeneity in dependent instruments (HEIDI) test.

Results

Characterization of bulk RNA-sequencing and functional enrichments

Following the methodology used to isolate uninfected (control) versus S. aureus-infected BMDMs, 3,775 DEGs were identified in the experimental group, comprising 1,686 upregulated and 2,089 downregulated genes. GO analysis revealed a network strongly enriched for innate immune activation, including positive regulation (GO:0045089) and activation (GO:0002218) of the innate immune response, reflecting can...

Discussion

By integrating TWAS and SMR analyses across translational datasets, this study identified a core network of genetic loci associated with osteomyelitis phenotypes. The identification of LYNX1 and MS4A14 as candidate gene targets represents a key contribution, providing mechanistic insight beyond conventional GWAS findings by incorporating transcriptomic context. This integrative approach highlights the potential of precision medicine strategies to advance the clinical management of severe bone infections.

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
clusterProfiler	Guangchuang Yu's Team (Sun Yat-sen University)	NA	Used for Gene Ontology (GO) functional enrichment analysis of differentially expressed genes.
DESeq2	Bioconductor Project (Open Source)	NA	Used to identify differentially expressed genes (DEGs)
GSEA / GseaVis	Broad Institute / Open Source Community	NA	Used for Gene Set Enrichment Analysis (GSEA) to identify coordinately expressed gene sets. GseaVis is likely used for visualization.
HEIDI (Heterogeneity in Dependent Instruments) test	Integrated function within the SMR software	NA	Used to test if the association observed in SMR is driven by linkage disequilibrium rather than a shared causal variant. P_HEIDI > 0.05 indicates no major heterogeneity
JTI (Joint-Tissue Imputation)	Zhou et al. (Open Source Model)	NA	A TWAS algorithm that estimates gene expression similarity and epigenetic chromatin accessibility to optimize prediction accuracy.
pheatmap	Raivo Kolde (R package)	NA	Used to generate heatmaps of differentially expressed genes for data visualization.
PrediXcan	Gamazon et al. (Open Source Model)	NA	A TWAS algorithm that applies elastic net regression with five-fold cross-validation to impute gene expression and test for association with disease risk.
SMR (Summary-data-based Mendelian Randomization)	Developed by a team at Fudan University (Open Source Software)	NA	Used to perform summary-based Mendelian randomization analysis to explore causal associations between gene expression and osteomyelitis risk. Version 1.0.3 was used.
STAR	Alexander Dobin (Cold Spring Harbor Laboratory)	NA	Used to align RNA-seq reads to the reference genome (mouse genome GRCm38/mm10).
UTMOST	Hu et al. (Open Source Model)	NA	A TWAS algorithm that enhances prediction accuracy by leveraging multi-tissue expression data using a sparse group LASSO framework.