This study employed a multi-stage analytical framework integrating genome-wide association studies, summary-data-based Mendelian randomization, tissue-specific expression quantitative trait loci, and Bayesian colocalization analyses. All analyses utilized publicly available summary-level GWAS data, eliminating the requirement for individual-level participant data and ethical approval.

Three independent gout GWAS datasets from European ancestry populations were included as the discovery cohort. The first dataset (Major TJ et al., n=642,075) was identified through the GWAS Catalog with identifier GCST90428600 (2). The second dataset (Carss K et al., n=458,440) from non-Finnish European individuals in the UK Biobank was designated GCST90474006 (10). The third dataset (Verma A et al., n=437,979) from the VA Million Veteran Program represented diverse disease-relevant endotypes with identifier GCST90475731 (11). To maximize statistical power and generate robust discovery associations, we conducted a fixed-effects inverse-variance weighted meta-analysis of these three GWAS datasets. The meta-analysis was performed using standard protocols, with heterogeneity assessed via Cochran’s Q test (p-value threshold >0.05 indicating homogeneity). For external validation, we utilized gout GWAS summary statistics from the FinnGen R12 consortium (n=327,457), a national biobank comprising well-characterized Finnish participants with comprehensive health registry data (12).

Circulating metabolite GWAS summary statistics were obtained from the metabolomic study by Karjalainen et al. (n=136,016 multi-ancestry participants), providing genome-wide associations for 233 metabolic biomarkers (13). The lipidomic GWAS dataset derived from Ottensmann et al. included 179 lipid species measured in 7,174 Finnish individuals from which 495 genetic associations were identified at novel and established loci (14). Plasma proteome GWAS data were obtained from two large-scale proteogenomic studies: the deCODE Genetics consortium proteome study (n=35,559 European participants) (15) and the UK Biobank Plasma Proteomic Project (UKB-PPP, n=54,306 European individuals), both providing cis-pQTL-disease associations (16).

We conducted two-sample MR analysis to evaluate the causal associations between circulating biomarkers and gout risk. The inverse-variance weighted (IVW) method was employed as the primary analytical approach, which assumes balanced horizontal pleiotropy and provides the most precise causal estimates when this assumption is satisfied. For exposures with more than one independent instrumental variable, we utilized either the Wald ratio method (for single SNP instruments) or the random-effect IVW method (for multiple SNPs). To evaluate the robustness of causal estimates and assess sensitivity to pleiotropy violations, we implemented several complementary MR methods including MR-Egger regression, and weighted median (WM) method. Results were considered robust when the direction and significance of causal effects remained consistent across multiple analytical methods.

To obtain comprehensive proteome-wide coverage while minimizing redundancy, we identified the intersection of proteins present in both the deCODE consortium proteomics dataset and the UK Biobank Pharma Protein Panel (UKB-PPP), resulting in 926 proteins subjected to MR analysis. This integrated approach leverages the strengths of both datasets and ensures robust protein-level associations with gout susceptibility.

Given the large number of exposures analyzed across three distinct omics categories—233 metabolites, 179 lipid species, and 926 proteins—we applied stringent multiple testing correction using the Benjamini-Hochberg false discovery rate (FDR) procedure within each omics category separately. This stratified FDR correction approach is particularly appropriate for multi-omics discovery studies, as it controls the false discovery rate while maintaining adequate statistical power for each biomarker class independently. For the discovery phase using the gout meta-GWAS, we considered associations statistically significant at an FDR-adjusted threshold of P_FDR < 0.05 within each omics category.

To establish robust and reliable associations, we employed a rigorous two-stage discovery and replication design. All biomarkers reaching significance in the discovery phase (gout meta-GWAS, FDR < 0.05) were systematically tested in an independent replication cohort (FinnGen R12 consortium gout GWAS). A finding was considered successfully replicated only when it satisfied both of the following criteria: (1) nominal statistical significance in the replication analysis (P < 0.05), and (2) consistent direction of effect between discovery and replication analyses (identical sign of effect estimate). Associations failing to meet both criteria were excluded from downstream functional analyses. This stringent validation approach ensures that identified causal biomarkers represent genuine, reproducible associations rather than chance findings.

Heterogeneity among SNP-specific causal estimates was assessed using Cochran’s Q test, with statistical significance indicating potential heterogeneity that might be attributable to pleiotropy or weak instrument bias. We evaluated horizontal pleiotropy—where genetic variants affect the outcome through pathways independent of the primary exposure—using MR-PRESSO (Mendelian Randomization Pleiotropy RESidual Sum and Outlier), which identifies and removes outlier variants that violate MR assumptions. The MR-Egger intercept test was additionally applied to estimate directional pleiotropy, where a statistically significant non-zero intercept (P < 0.05) indicates directional pleiotropy. Results from sensitivity analyses demonstrating substantial heterogeneity or pleiotropy were interpreted with caution and required additional validation through complementary methods.

For biomarkers successfully validated in replication analyses, we performed SMR analysis to identify specific genes whose expression is causally associated with gout risk, utilizing tissue-specific eQTL data from GTEx (kidney, liver, and whole blood). The heterogeneity in dependent instruments (HEIDI) test was applied as a critical quality control filter to distinguish genuine mediation effects from associations driven by linkage disequilibrium or shared variants between eQTL and GWAS summary statistics. We employed a stringent HEIDI test P-value threshold of P_HEIDI > 0.05 to exclude associations with evidence of linkage disequilibrium-driven false positives. Genes passing HEIDI filtering were retained for downstream colocalization analysis and functional validation.

To identify genomic loci where biomarker-associated variants, tissue-specific eQTL signals, and gout risk variants share a single causal variant, we performed multi-trait Bayesian colocalization analysis using the HyPrColoc (Hypothesis Prioritization for multi-trait Colocalization) framework. Unlike pairwise methods, this algorithm allows for the simultaneous assessment of colocalization across multiple traits by clustering them based on shared genetic drivers. We evaluated the posterior probability (PP) that the biomarker GWAS, tissue-specific eQTL, and gout meta-GWAS define a shared cluster of localized association. We set a stringent threshold of PP > 0.7 to indicate significant evidence of multi-trait colocalization. This integrative approach requires concordance at the variant level across all three molecular data types, substantially reducing the likelihood of spurious associations driven by linkage disequilibrium.

Candidate genes identified through both SMR analysis and colocalization were subjected to functional validation using quantitative reverse transcription PCR (qRT-PCR) and Western blot analysis in relevant tissue samples. For qRT-PCR, kidney tissue RNA was extracted using TRIzol reagent, and cDNA synthesis was performed using standard protocols. Real-time qRT-PCR was conducted in triplicate using PowerUp SYBR Green Master Mix, with thermal cycling optimized for each target. Relative gene expression was calculated using the 2^(-ΔΔCt) method normalized to housekeeping genes (GAPDH and ACTB). For Western blot analysis, total protein was extracted from kidney tissue lysates using RIPA buffer supplemented with protease and phosphatase inhibitors. Equal amounts of protein (30-50 μg) were resolved on polyacrylamide gels, transferred to PVDF membranes, and incubated with primary antibodies targeting candidate gene products. Protein bands were visualized using enhanced chemiluminescence and quantified by densitometry relative to β-actin or α-tubulin loading controls.

THP-1 cells (ATCC TIB-202) were maintained in RPMI 1640 medium (Gibco) supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin at 37 °C in 5% CO₂. For macrophage differentiation, THP-1 cells were seeded at 2 × 10⁵ cells/ml and treated with phorbol 12-myristate 13-acetate (PMA; 100 ng/ml) for 48 hours. Following differentiation, cells were rested in PMA-free medium for 24 hours prior to treatment. Differentiated THP-1 macrophages were then stimulated with monosodium urate (MSU) crystals at 200 μg/ml for 24 hours, with PBS-treated cells serving as vehicle controls. All experiments were conducted in triplicate with a minimum of three independent biological replicates.

To validate the transcriptional changes of the identified risk genes, total RNA was extracted from THP-1 macrophages treated with MSU crystals or vehicle control using TRIzol reagent (Invitrogen, USA) according to the manufacturer’s protocol. RNA concentration and purity were assessed using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific), with acceptable samples having A₂₆₀/₂₈₀ ratios between 1.8 and 2.0. Subsequently, 1 μg of total RNA was reverse-transcribed into cDNA using the PrimeScript RT Reagent Kit (Takara, Japan).

Quantitative real-time PCR was performed using SYBR Green Master Mix (Applied Biosystems) on an ABI 7500 Real-Time PCR System with the following thermal cycling conditions: initial denaturation at 95 °C for 10 minutes, followed by 40 cycles of 95 °C for 15 seconds and 60 °C for 60 seconds. Target genes assessed included five protein-coding genes—PRELID1, NIPAL1, LMAN2, and CAD—and one long non-coding RNA (lncRNA), AC093690.1. GAPDH and ACTB were employed as reference genes for normalization. Relative expression levels were calculated using the comparative 2^(–ΔΔCt) method and expressed as fold-change relative to vehicle-treated control cells. All primer sequences are listed in Table 1.

To verify the protein abundance of the candidate genes, total protein was extracted from THP-1 cells using RIPA lysis buffer (150 mM NaCl, 1% Triton X-100, 0.5% sodium deoxycholate, 0.1% SDS, 50 mM Tris-HCl pH 8.0) supplemented with protease and phosphatase inhibitor cocktails (Roche). Protein concentration was determined using the BCA Protein Assay Kit (Pierce Biotechnology).

Given the wide range of molecular weights among the target proteins (25 kDa to 240 kDa), appropriately-sized SDS-PAGE gels were selected for each target: 6% polyacrylamide gels for the high-molecular-weight protein CAD (240 kDa), and 10-12% gels for PRELID1 (25 kDa), NIPAL1 (34 kDa), and LMAN2 (35 kDa). Protein samples (30-50 μg per lane) were separated under denaturing conditions and transferred onto PVDF membranes (Millipore) using a wet transfer system.

Membranes were blocked with 5% non-fat milk in Tris-buffered saline with 0.1% Tween 20 (TBST) for 1 hour at room temperature and incubated overnight at 4 °C with primary antibodies against PRELID1 (1:1000; Proteintech, 10877-1-AP), NIPAL1 (1:1000; Sigma-Aldrich, HPA036765), LMAN2 (1:1000; Proteintech; 11496-1-AP), and CAD (1:1000; Proteintech; 16617-1-AP). β-actin (1:5000; Proteintech; 20536-1-AP) was used as the loading control. Following three 10-minute washes in TBST, membranes were incubated with HRP-conjugated secondary antibodies (1:5000) for 1 hour at room temperature. Protein bands were visualized using enhanced chemiluminescence (ECL) reagent and quantified using ImageJ software. Band intensities for each target protein were normalized to the corresponding β-actin band intensity on the same membrane, and results are expressed as fold-change relative to vehicle-treated control cells.

All statistical analyses were performed using R version 4.3.1 (R Foundation for Statistical Computing). Mendelian randomization analyses utilized the TwoSampleMR package (version 0.5.7) and Mendelian Randomization package (version 0.7.0). Instrumental variables were selected using genome-wide significance threshold (P < 5×10^-8), followed by LD clumping (r² < 0.001 within 10,000 kb window) using PLINK v1.90 with the 1000 Genomes Phase 3 European reference panel. Instrument strength was verified using F-statistic > 10. Summary-data-based Mendelian randomization was performed using SMR software version 1.3.1. HEIDI test was applied with default parameters to distinguish causal associations from linkage. Bayesian colocalization analysis was performed using the HyPrColoc package (version 1.0.0) with default prior probabilities. GWAS meta-analysis was conducted using METAL software with inverse-variance weighted fixed-effects model. Heterogeneity across studies was assessed using Cochran’s Q statistic and I² index. All visualization was performed using ggplot2 (version 3.4.4) and ComplexHeatmap (version 2.16.0) packages.

To establish a high-powered genetic foundation for identifying causal biomarkers, we performed a fixed-effects inverse-variance weighted meta-analysis of three large-scale gout GWAS datasets (N = 1,538,494). The final harmonized dataset covered 12,754,587 variants. We observed exceptional biological validity through the robust replication of established gout risk loci. Specifically, the classic missense variant in ABCG2 (rs2231142) and the lead variant in SLC2A9 (rs16890979) demonstrated associations exceeding genome-wide significance by orders of magnitude (P = 1.66×10⁻³⁰⁵ and P = 1.04×10⁻³⁰⁵, respectively), confirming the directionality and phenotype specificity of our meta-analysis. Furthermore, the genomic inflation factor (λGC) was 1.08, indicating that population stratification was effectively controlled and that the observed signals are driven by genuine polygenic architecture rather than systemic bias. We identified statistically significant heterogeneity in a subset of variants, which were rigorously excluded from downstream instrumental variable selection to ensure the robustness of the MR analyses. The resulting high-quality summary statistics served as the primary outcome for the subsequent multi-omics causal inference.

In MR analyses of 233 circulating metabolites, 32 showed significant causal associations with gout after FDR correction (P_FDR < 0.05), and all 32 replicated in FinnGen with consistent direction of effect (nominal P < 0.05). Most signals were lipid-related and clustered in triglyceride-rich very-low-density lipoprotein (VLDL) particles. Triglycerides and cholesterol components in medium and large VLDL were consistently associated with higher gout risk (e.g., triglycerides in medium VLDL: OR_discovery = 1.24, 95% CI 1.12–1.38; OR_validation = 1.19, 95% CI 1.05–1.34), and both VLDL particle concentrations and mean VLDL diameter were positively associated with gout. Similar patterns were observed for chylomicrons and extremely large VLDL, extending the lipoprotein–gout association across multiple particle classes (Figure 1).

Amino acid profiling identified circulating isoleucine as a strong causal risk factor, with larger effect sizes than lipid traits (OR_discovery = 1.86, 95% CI 1.25–2.77; OR_validation = 1.78, 95% CI 1.28–2.47), robust across multiple MR methods. In contrast, several lipoprotein compositional ratios showed protective associations, including the phospholipids-to-total lipids ratio in medium VLDL (OR_discovery = 0.86, 95% CI 0.78–0.96) and the free cholesterol-to-total lipids ratio in large HDL (OR_discovery = 0.84, 95% CI 0.73–0.96), suggesting that specific features of lipoprotein composition may mitigate gout risk (Figure 1).

Among 179 lipid species, multiple showed suggestive associations, but only triacylglycerol TAG (54:3) remained significant after FDR correction and successfully replicated in the independent cohort (OR_discovery = 1.13, 95% CI 1.07–1.19, FDR = 6.33×10^-5; OR_validation = 1.15, 95% CI 1.04–1.26). Sensitivity analyses and pleiotropy tests supported a causal effect. TAG (54:3), a marker of triglyceride-rich lipoproteins, is consistent with and likely reflects the broader VLDL-related metabolomic signals, while the overall lipidomic signal landscape remained sparse (Figure 2).

Proteome-wide MR of 926 plasma proteins identified two proteins with FDR-significant causal associations with gout, both protective: ISLR2 and ITIH3. Higher ISLR2 levels were associated with lower gout risk (OR_discovery = 0.89, 95% CI 0.83–0.95; OR_validation = 0.84, 95% CI 0.75–0.94), as were higher ITIH3 levels (OR_discovery = 0.86, 95% CI 0.79–0.94; OR_validation = 0.88, 95% CI 0.80–0.97). These associations were robust across alternative MR estimators, and pleiotropy/heterogeneity tests did not indicate major violations of MR assumptions, supporting a protective role of ISLR2 and ITIH3 in gout pathogenesis (Figure 3). Pleiotropy assessment for both proteins revealed minimal violations of MR assumptions.

To elucidate the transcriptomic architecture underlying the identified metabolic associations, we performed a multi-tissue SMR analysis, summarizing the prioritized effector genes and their tissue-specific expression patterns in Figure 4. The kidney cortex emerged as a key regulatory tissue, containing genes with distinct risk-increasing and protective profiles. PRELID1 (prelipin domain–containing 1) was identified as a significant renal risk gene that higher cortical expression was causally associated with increased gout risk (β_SMR = 0.050, P_SMR = 9.22×10^-5, P_HEIDI = 0.46) and with elevated levels of isoleucine and cholesteryl esters in large VLDL (β_SMR = 0.030, P_SMR = 2.97×10^-4), supporting a specific amino acid– and lipid-driven risk pathway. By contrast, we observed robust protective mechanisms in the kidney. CLEC18A (C-type lectin domain containing 18A) showed a significant inverse association with gout, whereby higher expression was linked to reduced disease susceptibility (β_SMR = −0.078, P_SMR = 1.95×10^-5, P_HEIDI = 0.50) and to lower circulating isoleucine levels. In addition, the long non-coding RNA RP11-392O17.1 displayed a protective profile that lower expression was causally associated with reduced isoleucine abundance (β_SMR = −0.032, P_SMR = 1.11×10^-4) and decreased gout risk (β_SMR = −0.034, P_SMR = 7.99×10^-4), implicating renal non-coding regulatory elements in metabolic pathogenesis.

Hepatic transcriptomic signals showed extensive colocalization with lipid-rich metabolic traits. NIPAL1 (NIPA-like domain containing 1) emerged as the top hepatic risk gene, with higher liver expression exerting a strong causal effect on gout risk (β_SMR = 0.074, P_SMR = 1.52×10^-8, P_HEIDI = 0.46) and on total lipids in chylomicrons and extremely large VLDL. The glucuronosyltransferase UGT2B17 provided particularly compelling evidence for a coordinated gene–metabolite–disease axis. Increased hepatic UGT2B17 expression was associated with both higher gout risk (β_SMR = 0.040, P_SMR = 3.96×^-5) and increased levels of serum total triglycerides and triglycerides in large VLDL (β_SMR = 0.038, P_SMR = 3.91×10^-6). This concordant directionality suggests that enhanced hepatic glucuronidation capacity may pathologically influence the abundance of triglyceride-rich lipoproteins, thereby predisposing individuals to gout.

In whole blood, we identified strong associations reflecting systemic metabolic regulation. The lncRNA AC093690.1 showed one of the largest effects in the study that higher expression was causally linked to a marked increase in gout risk (β_SMR = 0.155, P_SMR = 1.03×10^-9, P_HEIDI = 0.11) and was strongly associated with elevated serum isoleucine levels (β_SMR = 0.106, P_SMR = 1.77×10^-10). Additionally, LMAN2 (lectin, mannose-binding 2) was identified as a substantial risk factor (β_SMR = 0.339, P_SMR = 0.002), with its expression specifically linked to the mean diameter of VLDL particles, further supporting the systemic contribution of glycoprotein transport pathways to lipoprotein remodeling and disease etiology.

To definitively establish that the identified eQTLs, metabolic biomarkers, and gout risk are driven by shared genetic variants rather than coincidental overlap due to linkage disequilibrium, we performed multi-trait Bayesian colocalization analysis using the HyPrColoc framework. Strong evidence of colocalization was observed for our top prioritized candidate genes, reinforcing their roles as specific effectors of metabolic risk pathways. In the kidney cortex, we detected a robust colocalization signal at the PRELID1 locus (PP = 0.955). The variant rs6885410 was identified as the shared causal driver linking renal PRELID1 expression, alterations in phospholipid levels in small VLDL particles, and gout susceptibility, confirming PRELID1 as a high-confidence renal effector gene. In whole blood, LMAN2 exhibited extremely strong colocalization (PP = 0.945), with rs34604271 serving as the shared variant influencing the triglycerides-to-total-lipids ratio in very small VLDL. Similarly, CAD showed substantial evidence of colocalization (PP = 0.795) mediated by rs6547692, validating the mechanistic link between pyrimidine/purine metabolism, VLDL composition, and gout risk. Hepatic regulation was highlighted by NIPAL1, which demonstrated a high posterior probability of 0.886 (candidate SNP rs13146880) in the liver, linking hepatic expression to VLDL triglyceride ratios.

To verify the biological relevance of our computationally predicted effector genes, we assessed their transcriptional and translational responses in an in vitro model of acute gouty inflammation. Differentiated THP-1 macrophages were stimulated with MSU crystals (200 μg/mL), mimicking the pathogenic microenvironment of gout. qPCR analysis (Figure 5A) revealed significant transcriptional dysregulation among all five prioritized candidates upon MSU stimulation. Consistent with their identification as risk genes in our SMR analysis, the mRNA expression levels of PRELID1, NIPAL1, and LMAN2 were significantly upregulated in the MSU-treated group compared with vehicle controls (P < 0.05). Similarly, the long non-coding RNA AC093690.1 exhibited a marked increase in expression, validating its putative role as a pro-inflammatory regulator. Conversely, CAD, which was identified as a protective factor in our SMR analysis, displayed a significant downregulation in response to MSU crystals, further supporting its potential role in mitigating gout susceptibility.

We next sought to confirm whether these transcriptional changes translated to alterations in protein abundance. Western blot analysis of THP-1 cell lysates demonstrated concordance with the mRNA findings (Figure 5B). The protein levels of PRELID1, NIPAL1, and LMAN2 were substantially elevated following MSU exposure. In parallel, CAD protein abundance was significantly reduced in the inflammatory state. These experimental data provide crucial biological validation of our multi-omics framework, demonstrating that these genetically causal genes are dynamically regulated during the acute inflammatory response characteristic of gout pathogenesis.