We yielded summary statistics of 486 human blood metabolites from the metabolomics server (Shin et al., 2014), which was one of the most comprehensive studies of metabolites thus far. The association analysis was carried out for ∼2.1 million SNPs up on 7,824 individuals of European descent. After quality control, a total of 486 metabolites (i.e., 309 known and 177 unknown metabolites) were analyzed (Shin et al., 2014). We also obtained summary statistics of six CVDs also generated from individuals of European ancestry (Table 1), including AF (Nielsen et al., 2018), hypertension (measured via diastolic blood pressure [DBP], systolic blood pressure [SBP] and pulse pressure [PP] on patients with hypertension) (Evangelou et al., 2018), MI (Nikpay et al., 2015), CAD (Nikpay et al., 2015), AIS (Malik et al., 2018) and HF (Shah et al., 2019). Note that, we did not discover the evidence of overlapping subjects for GWASs of the six diseases. For all summary datasets, we performed the following quality control procedure: 1) deleted SNPs of non-biallelic; 2) excluded SNPs with no rs labels and duplicate SNPs; 3) excluded SNPs in the major histocompatibility complex region (chr6: 25.5–33.5 Mb); 4) retained SNPs not included in the 1000 Genome Project (The 1000 Genomes Project Consortium, 2015); 5) kept SNPs with minor allele frequency (MAF) > 0.01.

The Manhattan and QQ plots of p values for these diseases are shown in Supplementary Figures S1, S2, where an evident inflation in test statistics is observed. However, the estimated genomic control factor and the intercept of LDSC suggest the observed inflation is primarily due to polygenic signals rather than confounding influence such as population stratification and unknown cryptic relatedness (Supplementary Table S1). Therefore, we did not conduct genomic control on test statistics of any diseases and still employed the original summary datasets in our analysis.

For each metabolite we carefully selected a set of independent associated SNPs serving as candidate instrumental variables. To this aim, we applied the clumping procedure in PLINK (version v1.90b3.38) (Purcell et al., 2007). Following prior studies (Choi et al., 2019; Sanna et al., 2019; Yang et al., 2020), we set the primary and secondary significance levels of the index SNP at 1 × 10^–5, r ² to 0.1 and a physical distance of 500Kb, with the 1000 Genome Project as a reference panel. Due to the small sample size of metabolites, we here used a relatively relaxed statistical threshold of 1 × 10^–5 rather than the more stringent genome-wide significance level of 5 × 10^–8 (Sanna et al., 2019). Lower threshold would lead to few instrumental variables reserved for most of the analyzed metabolites; in contrast, higher threshold (e.g., 1 × 10^–5) is generally employed to generate more instrumental variables; therefore, larger variation of the exposure is explained, which has the potential to improve power in MR studies. Furthermore, to avoid the influence of horizontal pleiotropy, we relied on a conservative strategy by excluding some candidate instrumental variables that were located less than 1 Mb away from related loci of cardiovascular diseases and whose p values were less than 0.05 after the Bonferroni correction (Zeng and Zhou, 2019a; Zhao and Schooling, 2019). Intuitively, if a metabolite-associated SNP instrumental variable is also related to CVDs, then this instrument would be potentially invalid. Therefore, excluding such instrument would minimize the influence of horizontal pleiotropy.

Depending on selected instrumental variables of metabolites, we primarily applied the inverse-variance weighted (IVW) method to estimate their causal effects on cardiovascular diseases (Burgess et al., 2017; Bandres-Ciga et al., 2019). We deemed there was a statistically significant association if the estimated causal effect of a given metabolite had a false discovery rate (FDR) < 0.05. To assess the robustness of our results, we also performed several complimentary and sensitivity analyses: 1) the maximum likelihood method (Nguyen et al., 2015) as well as the weighted median-based method when instrumental variables might be invalid (Bowden et al., 2016); 2) the MR-Egger regression to evaluate the directional pleiotropy of instruments (Bowden et al., 2015); 3) the MR-PRESSO test to identify outliers (Verbanck et al., 2019); 4) the multivariable MR analysis to evaluate multiple metabolites simultaneously showing the association with the disease of focus (Burgess et al., 2013; Burgess and Thompson, 2017) (Supplementary Text); 5) the IVW analysis with the disease as an exposure and the metabolite as an outcome to examine the presence of reverse causality if a causal influence of the metabolite on the disease was identified. And these instrumental variables for CVDs were selected by a similar PLINK clumping procedure described above, but at a genome-wide significance level of 5 × 10^–8 (Supplementary Text).

To investigate whether the identified causal association between cardiovascular diseases and the metabolite can be attributable to common genetic foundation, we conducted the colocalization analysis using the R coloc package. We first extracted summary statistic information of SNPs within 50 Kb of an instrumental variable, and performed the colocalization analysis with default parameters. We then relied on only p values and minor allele frequencies to calculate five posterior probabilities (i.e., PP0, PP1, PP2, PP3 and PP4) (Giambartolomei et al., 2014). Among these, large PP3 indicates that both the disease and the metabolite are associated, but with different causal variants; while large PP4 (>80%) supports both the disease and the metabolite are associated and share a single causal variant (Giambartolomei et al., 2014; Steinberg et al., 2021).

For a significant causal association, we also conducted LDSC to study the genetic correlation between the disease and the metabolite with genome-wide SNPs (Bulik-Sullivan et al., 2015). Genetic correlation provides an overall perspective into shared genetic foundation underlying the two types of phenotypes (van Rheenen et al., 2019).

Finally, based on all metabolites showing suggestively significant association with any of the six cardiovascular diseases, we performed a metabolic pathway analysis using MetaboAnalyst5.0 (Chong et al., 2018). The metabolic pathway analysis includes two datasets: the Small Molecular Pathway database (SMPDB) (Frolkis et al., 2010) and the KEGG database (Kanehisa et al., 2012).

The number of instrumental variables for metabolites ranged from 3 to 631, with a median number of 22. On average, the selected SNPs explained 10.1% of phenotypic variance across all the 486 metabolites (Supplementary Figure S3). Importantly, the minimum F statistics were above 10 (ranging from 17.4 to 24.9) (Supplementary Table S2), indicating that weak instrumental bias is unlikely to occur (Burgess et al., 2017). Using these instrumental variables, we assessed the causal association between 486 metabolites and six cardiovascular diseases, and identified a total of 310 suggestive associations (p < 0.05; corresponding to 207 unique metabolites), including 198 associations for 135 known metabolites and 112 associations for 72 unknown metabolites (Supplementary Table S3). Among these, there were 20, 25, 30, 26, 26, 25, 24, and 22 associations known metabolites (Figure 1) and 19, 11, 18, 12, 16, 18, 6, and 12 associations unknown metabolites (Supplementary Figure S4) related to AF, CAD, DBP, PP, SBP, HF, AIS, and MI, respectively. After the multiple-testing correction, we obtained four significant associations (FDR < 0.05): bradykinin, des-arg(9) (odds ratio [OR] = 1.160, 95% confidence intervals [CIs]: 1.080–1.246, FDR = 0.022) on AIS, N-acetylglycine (OR = 0.946, 95%CIs: 0.920–0.973, FDR = 0.023) on DBP, X-09026 (OR = 0.845, 95%CIs: 0.779–0.916, FDR = 0.021) on DBP and X-14473 (OR = 0. 938, 95%CIs = 0.907–0.971, FDR = 0.040) on DBP (Supplementary Table S3). In addition, there are 14 promising associations (0.05 < FDR < 0.10), such as tryptophan betaine (OR = 0.882, 95%CIs: 0.827–0.940, FDR = 0.058) on AF and N-acetylornithine (OR = 0.860, 95%CIs: 0.797–0.929, FDR = 0.055) on CAD (Supplementary Table S3).

We observe that bradykinin, des-arg (9) also shows a suggestive association with DBP (OR = 1.020, 95%CIs: 1.007–1.032, p = 0.002) and PP (OR = 0.938, 95%CIs: 0.907–0.971, p = 0.034), and that N-acetylglycine (OR = 0.960, 95%CI: 0.936–0.985, p = 0.002) and X-14473 (OR = 0.967, 95%CI: 0.943–0.992, p = 0.010) are suggestively associated with DBP (Supplementary Table S4). Moreover, 52 metabolites are associated with at least two cardiovascular diseases (p < 0.05). Interestingly, some of these metabolites had the opposite causal effect across the diseases, such as a positive influence of linolenate [alpha or gamma; (18:3n3 or 6)] on HF (OR = 1.566) but a negative impact on AIS (OR = 0.648) (Supplementary Table S4), implying distinct functional roles of metabolites in the development of cardiovascular diseases.

The full results of sensitivity analyses are shown in Supplementary Table S5. Generally, the weighted median method and the maximum likelihood method generate similar causal effect estimates compared to the fixed-effects IVW MR method. The results of sensitivity analyses for the four significant metabolites on AIS/DBP are summarized in Figure 2. Again, the two methods show robust causal associations, such as bradykinin, des-arg(9) on AIS (P _{Weight-median} = 0.0003 and P _Likelihood = 0.0001), N-acetylglycine on DBP(P _{Weight-median} = 0.0027 and P _Likelihood = 0.0002), X-09026 on DBP (P _{Weight-median} = 0.0017 and P _Likelihood = 0.0012) and X-14473 on DBP (P _{Weight-median} = 0.0006 and P _Likelihood = 0.0015). The intercept of MR-Egger regression does not deviate significantly from zero, indicating the absence of horizontal pleiotropy; however, the causal effects of N-acetylglycine (P _MR-Egger = 0.267), X-09026 (P _MR-Egger = 0.770) and X-14473 (P _MR-Egger = 0.394) on DBP were nonsignificant in terms of the MR-Egger test, in line with the prior finding that the MR-Egger method is in general less efficient compared to other used methods. As a further sensitivity analysis, for the four significant metabolites we tried a more stringent significance level of 5 × 10^–8 to screen instrumental variables. Unfortunately, only bradykinin, des-arg(9) and N-acetylglycine had SNP instruments at this level. Nevertheless, the corresponding associations are still significant, with P _{Weight-median} = 0.0007 and P _Likelihood = 0.0070 for bradykinin, des-arg(9) on AIS and P _{Weight-median} = 8.23 × 10^–10 and P _Likelihood = 5.41 × 10^–9 for N-acetylglycine on DBP, respectively. Moreover, using SNP instruments obtained at the level of 1 × 10^–6, we also produce a significant association between X-09026 and DBP (P _{Weight-median} = 0.0013 and P _Likelihood = 0.0004). It is easy to see that all of these results are highly consistent with those obtained using instruments with a relatively relaxed significance level of 1 × 10^–5 that we applied in our main analysis, indicating the robustness of these identified association signals. In addition, we created scatter plots of SNP effect sizes for the four metabolites and AIS/DBP, and show that no instrumental variables behave as potential outliers. Funnel plots based on individual causal effect estimates of metabolites on AIS/DBP display a symmetrical pattern and provide little evidence of horizontal pleiotropy (Supplementary Figures S5–S8). MR-PRESSO also does not support the presence of horizontal pleiotropy and instrumental outliers (P _outlier > 0.05).

We applied the multivariate MR method to analyze whether the causal effect of one metabolite on DBP would be affected by other metabolites. It is shown that the causal effects estimated with the multivariate MR method are consistent with the unadjusted ones obtained via the fixed-effects IVW MR method for three identified metabolites (i.e., N-acetylglycine, X-09026 and X-14473) (Supplementary Table S6), implying the independent role of these metabolites. We also examined the causal relationship between the four identified metabolites and observe a significant causal effect of X-09026 on X-14473 (β = 0.101, 95%CI: 0.041–0.160, p = 9.03 × 10^–4) and N-acetylglycine on X-14473 (β = -0.119, 95%CI: 0.235 ∼ -0.003, p = 0.045) (Supplementary Table S7). Interestingly, we find that X-09026 and n-acetylglycine exhibit a causal effect on X-14473; however, these three metabolites seem independently affect DBP (Figure 3).

We next carried out a reverse causality analysis by using instrumental variables of AIS/DBP to make causal inference for metabolites. To this aim, using the similar clumping procedure in PLINK we selected 346 and 28 independent index SNPs (p < 5 × 10^–8) as instrumental variables of DBP or AIS and carried out the fixed-effects IVW MR estimation. But we find little evidence supporting the presence of reverse causal relationship between these four metabolites and AIS/DBP (Supplementary Table S8).

To identify whether the association between metabolites and cardiovascular diseases was attributable to common genetic component, we conducted the LDSC and colocalization analyses. The LDSC analysis shows non-significant genetic correlations between AIS and bradykinin, des-arg(9) (r _g = 0.104, p = 0.662), DBP and N-acetylglycine (r _g = −0.043, p = 0.615), and DBP and X-14473 (r _g = 0.230, p = 0.468) (Supplementary Table S9). We here highlight that the genetic correlation between DBP and X-09026 cannot be calculated because the heritability of the metabolite X-09026 is negatively estimated. In addition, we evaluated whether genome-wide SNPs of four metabolites that were significantly associated with AIS/DBP were co-localized. The colocalization results suggest that none of previously detected signals is responsible for a single shared genetic variant (PP4 < 80%) (Supplementary Table S10), suggesting the identified associations are unlikely due to confounding factors by LD or common causal SNPs.

We identified a total of 15 important metabolic pathways involved in the pathogenesis of six cardiovascular diseases by all the identified metabolites (Table 2). The results show that some diseases have common metabolic pathways, such as “Caffeine metabolism” shared by CAD (p = 0.002), SBP (p = 0.033) and MI (p = 0.044), “Aminoacyl-tRNA biosynthesis” shared by AD (p = 0.003) and HF (p = 0.004), “Carnitine Synthesis” shared by AF (p = 0.019) and CAD (p = 0.019), “Bile Acid Biosynthesis” shared by CAD (p = 0.044) and HF (p = 0.044), “Primary bile acid biosynthesis” shared by DBP (p = 0.027) and HF (p = 0.040). In addition, we also find that both Glycine and l-cysteine in HF underwent three pathways, which are “Aminoacyl-tRNA biosynthesis” (p = 0.004), “Glutathione metabolism” (p = 0.016) and “Glycine, serine and threonine metabolism” (p = 0.021). These findings offer further insight into the metabolic mechanism for cardiovascular diseases.