In this study, we selected protein quantitative trait loci (pQTL) data from two large-scale plasma proteomics studies (deCODE Genetics study and the UK Biobank Pharma Proteomics Project (UKB-PPP)) as exposure variables and used genome-wide association study (GWAS) data from two cohorts (UK Biobank and the FinnGen cohort) as outcome measures to investigate the potential causal relationship between plasma proteins and AML risk. To avoid bias caused by population overlap (UK biobank and UKB-PPP), we conducted three independent MR analyses (deCODE to FinnGen, deCODE to UK Biobank, and UKB-PPP to FinnGen). We identified proteins that were nominally significant in at least two analyses and aligned in direction in three analyses as potential causal biomarkers. To assess the robustness of our findings, Bayesian colocalization analysis was performed. To further examine the expression patterns of the candidate biomarker genes across different cell types and reveal their cell-specific roles in AML, we utilized single-cell RNA sequencing (scRNA). Meanwhile, we analyzed the expression patterns of the candidate biomarker genes in AML subtypes defined by ELN and WHO guidelines (17, 18). The candidate proteins were then experimentally validated using enzyme-linked immunosorbent assay (ELISA). Additionally, a phenome-wide association study (PheWAS) was conducted to confirm the safety of the identified targets. Finally, we carried out druggability assessments to evaluate the potential of the identified plasma proteins as therapeutic targets. The detailed study design workflow is illustrated in Figure 1 .

GWAS data for AML were obtained from two datasets: the UK Biobank (cases = 220, controls = 456,128) and the FinnGen cohort (cases = 244, controls = 314,192). pQTL data were sourced from two studies: the Icelandic deCODE Genetics study (deCODE) and the UK Biobank Pharma Proteomics Project (UKB-PPP). The deCODE study provided pQTL data for 4,907 plasma proteins, measured in 35,559 Icelandic individuals using the SomaScan v4 platform. The UKB-PPP dataset included detailed pQTL mapping for 4,572 proteins, measured in 54,219 participants from the UK Biobank using the Olink platform.

To identify pQTLs, we applied the following criteria: To identify pQTLs, we applied the following criteria (19): (1) The SNP was located within a vicinity of ±1 Mb around the gene region (cis-acting pQTLs); (2) The SNP-protein association reached a genome-wide significant threshold of P < 5 × 10^-8; (3) The independence assumption was satisfied through linkage disequilibrium clumping (r² < 0.001); (4) Palindromic SNPs were excluded to avoid allele orientation ambiguity; (5) The SNPs with an F-statistic < 10 were excluded to ensure strong instrument strength and minimize weak instrument bias. Detailed results for each protein, including the number of instruments, variance explained (R²), and F-statistic are provided in the Supplementary Table 1 .

In this study, we conducted MR analyses using plasma proteins as exposure variables and AML as the outcome variable. The selection criteria for pQTLs strictly adhered to the standards described earlier. For proteins with three or more SNPs, the inverse-variance weighted (IVW) method was used as the primary approach. In addition, sensitivity analyses, including MR-Egger, weighted median, and weighted mode were conducted to evaluate pleiotropy and heterogeneity. For proteins instrumented by a single SNP, the Wald ratio method was applied. And for proteins instrumented by two or less SNP, sensitivity analyses were not performed (20). We further applied Steiger filtering to verify the direction of causality (exposure to outcome). Furthermore, to rule out potential reverse causation, we performed reverse MR analysis for all significant proteins identified in the forward analysis, using AML as the exposure and the respective protein as the outcome. For the reverse MR analysis, instrumental variables for AML were selected using a genome-wide significance threshold of P<5×10 ^-8. If no sufficient instruments were available at this threshold, a more lenient threshold of P<5×10^–6 was applied to obtain an adequate number of SNPs for analysis. To account for multiple testing, the false discovery rate (FDR) correction was applied; a causal association was considered significant if the adjusted p-value (PFDR) < 0.05, and suggestive evidence was noted for nominal significance (P < 0.05). All analyses were performed using the R package “TwoSampleMR” (v0.5.6).

We performed enrichment analysis on the plasma proteins identified through MR to investigate their potential biological functions and involvement in pathways. All analyses were conducted in metascape online tool (https://metascape.org) and visualizing the results.

Bayesian colocalization analysis was performed using the “coloc” package with default parameters to estimate the probability that two traits share the same causal variant. This approach helps determine whether observed associations are due to causal effects of genetic variants on traits, rather than being influenced by LD or other confounding factors. The colocalization analysis evaluates five hypotheses: H0: Neither the exposure nor the outcome is associated with the genomic region; H1: The exposure is associated with the genomic region, but the outcome is not; H2: The outcome is associated with the genomic region, but the exposure is not; H3: The exposure and the outcome are associated with the genomic region through different SNPs; H4: The exposure and the outcome are associated with the genomic region through a shared SNP. We calculated the posterior probabilities (PP) for each hypothesis and determined the presence of colocalization evidence for a protein based on the condition PPH3 + PPH4 > 0.5 (21).

Firstly, we utilized UniPort (https://www.uniprot.org/) to investigate the structure and cellular localization of proteins, and analyzed the tissue-specific expression of proteins through GeneCards (https://www.genecards.org/). Then, we analyzed the scRNA data through the Single Cell Portal database (https://singlecell.broadinstitute.org/) to evaluate the gene expression differences of TCL1A, TNFSF8, WFDC1 and TNFAIP8 in different cell types. The access number of the dataset is SCP1987, which includes bone marrow mononuclear cells from 42 patients with acute myeloid leukemia and 10 healthy individuals without specific immune cell enrichment or removal of blast cells. This dataset contains 254,910 cells and 33,947 genes (22). Finally, we analyzed the gene expression patterns of each subtype of AML based on the ELN and WHO guidelines using the BeatAML database (http://www.vizome.org/aml2/). Group comparisons were performed using one-way ANOVA or Kruskal-Wallis tests based on data distribution normality, with appropriate post-hoc tests (Tukey’s or Dunn’s test) for multiple comparisons.

To validate the robustness of the MR results, we collected serum samples from 6 age and sex-matched healthy controls and 20 AML patients (12 newly diagnosed patients, 4 in complete remission (CR) after treatment, and 4 with relapse following treatment, the detail information was list in Supplementary Table 2 ). The TNFSF8 ELISA kit (catalog number EH0120) was purchased from FineTest (https://www.fn-test.com/, China), the TNFAIP8 kit (catalog number ELH-TNFAIP8) was obtained from RayBiotech (https://www.raybiotech.com/, USA); and the WFDC1 and TCL1A kits (catalog numbers E9351h and E1709h, respectively) were sourced from EiAab (https://www.eiaab.com.cn/, China). The levels of TCL1A, TNFSF8, WFDC1, and TNFAIP8 in the serum were measured according to the manufacturer’s instructions for each respective kit. Briefly, all serum sample measurements were performed in triplicate. The mean value was used for subsequent analysis unless the coefficient of variation (CV) exceeded 15%, in which case the measurement was repeated. The intra-assay and inter-assay CVs for each target protein were below 10% and 15%, respectively. Standard curves were generated for each assay, all of which demonstrated excellent linearity (R² > 0.99). The limits of detection and quantification were determined for each kit according to the manufacturer’s protocols. All laboratory personnel were blinded to the clinical status of the samples during the measurement process. Any outliers were identified using the ROUT method (Q = 1%) and were excluded from the final analysis. Group comparisons were performed using one-way ANOVA or Kruskal-Wallis tests based on data distribution normality, with appropriate post-hoc tests (e.g., Tukey’s or Dunn’s test) for multiple comparisons.

Clinical data from 12 newly diagnosed AML patients were collected to analyze the correlation between plasma protein levels of TNFAIP8, TCL1A, WFDC1, and TNFSF8 and the clinical characteristics of the patients. Pearson correlation analysis was performed, and a p-value < 0.05 was considered statistically significant.

PheWAS was performed to evaluate the pleiotropic effects of potential therapeutic targets and possible adverse effects. The outcome involved obtaining phenotypic data from the Finnish database in version R10, encompassing 2408 phenotypes categorized into 46 groups. This extensive dataset was employed for phenome-wide MR analysis. P < 0.05 was considered statistically significant.

The target genes were submitted to the Drug Signatures Database (DSigDB, http://dsigdb.tanlab.org/DSigDBv1.0/) to explore potential interactions between the proteins and available drugs. The target genes were input into the Enrichr gene set enrichment analysis platform (https://maayanlab.cloud/modEnrichr/) to access the DSigDB database and predict potential drug candidates that could interact with the genes of interest.

To systematically identify therapeutic targets, we integrated multi-center data from exposure-outcome pairs. Three independent analyses identified 65, 69, and 86 AML-associated plasma proteins (p<0.05) from deCODE to FinnGen, deCODE to UK Biobank, and UKB-PPP to FinnGen, respectively ( Supplementary Figures S1A-C ; Supplementary Table 3 ). Enrichment analysis revealed that these proteins were associated with immune and inflammatory processes ( Supplementary Figures S1D-F ). We then screened for proteins consistently detected in at least two analyses with the same directional effects across all three datasets. This approach identified six candidate proteins: GRAP2, TNFAIP8, TCL1A, WFDC1, C7, and TNFSF8 ( Figure 2 ). In addition, in the reverse MR analysis for six candidate proteins, no causal association were found ( Supplementary Table 4 ). Notably, five proteins exhibited risk-promoting effects for AML progression (GRAP2 OR = 9.32, 95% CI: 2.28-38.151; TNFAIP8 OR = 3.66, 95% CI: 1.69-7.90; TCL1A OR = 2.34, 95% CI: 1.44-3.78; WFDC1 OR = 1.85, 95% CI: 1.40-2.46; C7 OR = 1.29, 95% CI: 1.10-1.52);, while TNFSF8 showed protective potential (OR = 0.33, 95% CI: 0.17-0.67). Further, we compared the top cis-pQTLs (the most significant SNPs associated with protein levels) for each key protein between the deCODE (SomaScan) and UKB-PPP (Olink) datasets. We observed high consistency in the top cis-pQTLs for four of the six proteins: WFDC1 (rs400345 in both platforms), TCL1A (rs78986913 in both platforms), TNFAIP8 (rs1035376 in both platforms) GRAP2 (rs148328786 in both platforms). For TNFSF8 and C7, the top cis-pQTL differed between platforms (TNFSF8, rs1006026 in deCODE vs. rs10081728 in UKB-PPP) (C7, rs72758315 in UKB-PPP, rs79534924 in deCODE).

To validate the genetic causality of these candidate proteins, we performed Bayesian colocalization analysis. The results suggested that TCL1A (PPH3+PPH4 = 59.65, 51.22, 64.99), TNFSF8 (PPH3+PPH4 = 68.62, 61.06, 71.63), WFDC1 (PPH3+PPH4 = 52.10, 65.62, 62.41), and TNFAIP8 (PPH3+PPH4 = 55.46, 53.92, 55.35) shared causal variants with AML across all three analyses (deCODE to FinnGen, deCODE to UK Biobank, and UKB-PPP to FinnGen, Figure 3 ). However, the PPH3+PPH4 values for GRAP2 (49.49, 38.80, 38.39) and C7 (53.06, 46.86, 51.32) indicated weak colocalization with AML in the genome. The detail results of colocalization analysis were list in Supplementary Table 5 .

Having established genetic associations, we next investigated the biological relevance of these proteins. The structure of protein was shown in Supplementary Figure S2A , and subcellular localization of the protein indicated that TNFSF8 is a cell-surface receptor (membrane), WFDC1 is a secreted protein, TCL1A is an intracellular protein (cytoplasm; endoplasmic reticulum; microsome; nucleus), and TNFAIP8 is also an intracellular protein (cytoplasm) ( Supplementary Figure S2B ). Tissue expression profiles showed that WFDC1, TCL1A and TNFAIP8 were mainly localized in immune cells and blood components, unfortunately, there is no data on TNFSF8 protein in GeneCards ( Supplementary Figure S3 ). Moreover, scRNA-seq results revealed that WFDC1 was mainly detected in GMP, MEP and MPP; TCL1A was mainly expressed in Pre-B cells and Pro-B cells; TNFAIP8 was primarily expressed in CD4+ T cells, CD8+ T cells, while TNFSF8 was predominantly expressed in CD4+ T cells ( Figure 4A ). Finally, we analyzed the expression of these genes in the AML subtypes defined by the ELN and WHO guidelines. The expression of TNFAIP8 was related to ELN2017 ( Figure 4B ). Their expressions varied among different FAB subtypes ( Figure 4C ). WFDC1 might be specific to RUNX1-RUNXT1, and TNFSF8 might be specific to CBFB-MYH11 ( Figure 4D ). TCL1A had lower expression in FLT3-ITD-positive patients, WFDC1 and TNFSF8 had lower expression in NPM1-positive patients, while TNAFIP8 had higher expression in both FLT3-ITD and NPM1-positive patients ( Figures 4E, F ).

We further verified the above results using ELISA, and ELISA results confirmed significantly elevated expression levels of TNFAIP8, WFDC1 and TCL1A, and low TNFSF8 expression in AML patients compared to healthy controls (P<0.05). Notably, in AML patients who achieved CR after treatment, the expression levels of TNFAIP8 and WFDC1 decreased, the expression levels of TNFSF8 increased, while TCL1A showed no significant change. Furthermore, TNFAIP8 levels were significantly higher and TNFSF8 levels were significantly lower in relapsed patients compared to CR patients, while the other proteins showed no significant differences ( Figure 4G ). Finally, to assess clinical significance of these proteins, we collected clinical data from 12 newly diagnosed AML patients and analyzed the correlation between plasma TNFAIP8, TCL1A, WFDC1, and TNFSF8 protein levels and patient clinical features. However, no significant correlations were found between the expression levels of these plasma proteins and the clinical parameters ( Supplementary Figure S4 ).

The above results indicated these four proteins could use as therapeutic target; to assess the potential beneficial or harmful effects of the four AML-associated plasma proteins on other phenotypes, we performed a phenome-wide association study (PheWAS). The results were shown in Supplementary Figure S5 and Supplementary Table 6 . These associations highlight both therapeutic potential and possible off-target effects of modulating these proteins.

Furthermore, to evaluate the potential of the identified proteins as drug targets, we investigated their interactions with known drugs using the DSigDB database. We found 13 drugs associated with TNFAIP8, 6 with TNFSF8, and 4 with WFDC1. However, no related drugs were identified for TCL1A ( Figure 5 ; Supplementary Table 6 ).