A Mendelian randomization study of the role of lipoprotein subfractions in coronary artery disease

Table 1 describes all GWAS summary datasets used in this study, including two GWAS of the traditional lipid risk factors (Willer et al., 2013; Hoffmann et al., 2018), two recent GWAS of the human lipidome (Kettunen et al., 2016; Davis et al., 2017), and three GWAS of CAD or MI (Nikpay et al., 2015; Nelson et al., 2017; Abbott et al., 2018). In the two GWAS of the lipidome (Kettunen et al., 2016; Davis et al., 2017), high-throughput nuclear magnetic resonance (NMR) spectroscopy was used to measure the circulating lipid and lipoprotein traits (Soininen et al., 2009). We investigated the 82 lipid and lipoprotein traits measured in these studies that are related to very-low-density lipoprotein (VLDL), LDL, intermediate-density lipoprotein (IDL), and HDL subfractions and particle sizes. All the subfraction traits are named with three components that are separated by hyphens: the first component indicates the size (XS, S, M, L, XL, XXL); the second component indicates the fraction according to the lipoprotein density (VLDL, LDL, IDL, HDL); the third component indicates the measurement (C for total cholesterol, CE for cholesterol esters, FC for free cholesterol, L for total lipids, P for particle concentration, PL for phospholipids, TG for triglycerides). For example, M-HDL-P refers to the concentration of medium HDL particles.

Aside from the concentration and content of lipoprotein subfractions, the two lipidome GWAS also measured the traditional lipid traits (TG, LDL-C, HDL-C), the average diameter of the fractions (VLDL-D, LDL-D, HDL-D) and the concentration of apolipoprotein A1 (ApoA1) and apolipoprotein B (ApoB). A full list of the lipoprotein measurements investigated in this article can be found in Appendix 1.

Genetic correlation is a measure of association between the genetic determinants of two phenotypes. It is conceptually different from epidemiological correlation that can be directly estimated from cross-sectional data. In this study, we applied the LD-score regression (Bulik-Sullivan et al., 2015) to the lipidome GWAS (Kettunen et al., 2016; Davis et al., 2017) to estimate the genetic correlations between the lipoprotein subfractions, particle sizes, and traditional risk factors. We then removed lipoprotein subfractions and particle sizes that are strongly correlated with the traditional risk factors, defined as an estimated genetic correlation > 0.8 with TG, LDL-C, HDL-C, ApoB, or ApoA1 in the GWAS published by Davis et al., 2017. Because these traits are largely co-determined with the traditional risk factors, they do not represent independent biological mechanisms and may lead to multicollinearity issues in multivariate MR analyses. Finally, we obtained an independent estimate of the genetic correlations between the selected traits by applying the LD score regression to the GWAS published by Kettunen et al., 2016. We used Bonferroni's procedure to correct for multiple testing (familywise error rate at 0.05).

For MR, we employed a three-sample design (Zhao et al., 2019b) in which one GWAS was used to select independent genetic instruments that are associated with one or several lipoprotein measures. The other two GWAS were then used to obtain summary associations of the selected SNPs with the exposure and the outcome, as in a typical two-sample MR design (Pierce and Burgess, 2013; Hemani et al., 2016). More specifically, the selection GWAS was used to create a set of SNPs that are in linkage equilibrium with each other in a reference panel (distance >10 megabase pairs, $r^2<0.001$). This was done by ordering the SNPs by the p-values of their association with the trait(s) under investigation and then selecting them greedily using the linkage-disequilibrium (LD) clumping function in the PLINK software package (Purcell et al., 2007). To avoid winner's curse, we require the other two GWAS to have no overlapping sample with the selection GWAS.

As the GWAS published by Davis et al., 2017 has a smaller sample size, we used it to select the genetic instruments so the larger dataset can be used for statistical estimation. In univariable MR, associations of the selected SNPs with the exposure trait (a lipoprotein subfraction or a particle size trait) were obtained from the GWAS published by Kettunen et al., 2016 and the associations with MI were obtained using summary data from an interim release of UK BioBank (Abbott et al., 2018). To maximize the statistical power, we used the so-called ‘genome-wide MR’ design. Independent SNPs are selected by using LD clumping, but we do not truncate the list of SNPs by their p-values. More details about this design can be found in a previous methodological article (Zhao et al., 2019b).

To control for potential pleiotropic effects via the traditional risk factors, we performed two multivariable MR analyses for each lipoprotein subfraction or particle size under investigation. The first multivariable MR analysis considers four exposures: TG, LDL-C, HDL-C, and the lipoprotein measurement under investigation. The second multivariable MR analysis replaces LDL-C and HDL-C with ApoB and ApoA1, in accordance with some recent studies (Richardson et al., 2020). SNPs were ranked by their minimum p-values with the four exposures and are selected as instruments only if they were associated with at least one of the four exposures (p-value $\le {10}^{-4}$). Both multivariable MR analyses used the Davis (Davis et al., 2017) and GERA (Hoffmann et al., 2018) datasets for instrument selection, the Kettunen (Kettunen et al., 2016) and GLGC (Willer et al., 2013) datasets for the associations of the instruments with the exposures, and the CARDIoGRAMplusC4D + UK Biobank (Nelson et al., 2017) dataset for the associations with CAD.

For univariable MR, we used the robust adjusted profile score (RAPS) because it is more efficient and robust than many conventional methods (Zhao et al., 2020; Zhao et al., 2019b). RAPS can consistently estimate the causal effect even when some of the genetic variants violate instrumental variables assumptions. For multivariable MR, we used an extension to RAPS called GRAPPLE to obtain the causal effect estimates of multiple exposures (Wang et al., 2020). GRAPPLE also allows the exposure GWAS to have overlapping sample with the outcome GWAS, while the original RAPS does not. We assessed the strength of the instruments using the modified Cochran's Q statistic (Sanderson et al., 2019). Because many lipoprotein subfraction traits were analyzed simultaneously, we used the Benjamini-Hochberg procedure to correct for multiple testing (Benjamini and Hochberg, 1995) and the false discovery rate was set to be 0.05. More detail about the statistical methods can be found in Appendix 3.

To obtain genetic markers, we selected SNPs that are associated with the lipoprotein measurements identified in the MR (p-value $\le 5\times {10}^{-8}$) and CAD (p-value $\le 0.05$) but are not associated with LDL-C or ApoB (p-value $\ge {10}^{-3}$). To maximize the power of this exploratory analysis, we meta-analyzed the results of the two lipidome GWAS (Kettunen et al., 2016; Davis et al., 2017) by inverse-variance weighting. For the associations with LDL-C and CAD, we used the GWAS summary data reported by the GLGC (Willer et al., 2013) and CARDIoGRAMplusC4D (Nelson et al., 2017) consortia. We used LD clumping to obtain independent markers (Purcell et al., 2007) and then validate the markers using tissue-specific gene expression data from the GTEx project.

Because we had multiple GWAS summary datasets for the lipoprotein subfractions and CAD/MI (Table 1), we swapped the roles of the GWAS datasets in the three-sample MR design whenever permitted by the statistical methods to obtain multiple statistical estimates. These estimates are not completely independent of the primary results, but they can nonetheless be used to assess replicability. As a sensitivity analysis, We further analyzed univariable MR using inverse-variance weighting (IVW) (Burgess et al., 2013) and weighted median (Bowden et al., 2016) and compared with the primary results obtained by RAPS. We also assessed the assumptions made by RAPS using some diagnostic plots suggested in previous methodological articles (Zhao et al., 2019b).

We obtained the genetic correlations of the lipoprotein subfractions and particle sizes with the traditional lipid risk factors: TG, LDL-C, HDL-C, ApoB, and ApoA1 (Table 1). We found that almost all VLDL subfractions traits (besides those related to very small VLDL subfraction) and the mean VLDL particle diameter have an estimated genetic correlation with TG very close to 1. Most traits related to the large and very large HDL subfractions also have a high genetic correlation with HDL-C and ApoA1.

After removing traits that are strongly correlated with the traditional risk factors, we obtained 27 traits that may involve independent genetic mechanisms. Figure 1 shows the genetic correlation matrix for these traits and the traditional lipid factors. The selected traits can be divided into two groups based on whether they are related to VLDL/LDL/IDL particles or HDL particles. Within each group, most traits were strongly correlated with the others. In the first group, most traits had a positive genetic correlation with LDL-C and ApoB, while in the second group, most traits had a positive genetic correlation with HDL-C and ApoA1. Exceptions include LDL-D, which had a negative but statistically non-significant genetic correlation with LDL-C and ApoB, and S-HDL-P and S-HDL-L, which showed no or weak genetic correlation with HDL-C and ApoA1.

Figure 1

Download asset Open asset

Figure 2 shows the estimated causal effect of the selected lipoprotein measurements on MI or CAD that are statistically significant (false discovery rate = 0.05). The unfiltered results can be found in Appendix 3, which also contains results of the sensitivity and replicability analyses.

Figure 2

Download asset Open asset

The concentration and lipid content of VLDL, LDL, and IDL subfractions showed harmful and nearly uniform effects on MI in univariable MR. However, after adjusting for the traditional lipid risk factors, the effects of these ApoB-related subfractions become close to zero (besides IDL-FC in one multivariable analysis). The mean diameter of LDL particles (LDL-D) showed a harmful effect on MI in univariable MR, though the effect was smaller than those of the LDL subfractions in univariable MR. The estimated effect of LDL-D was attenuated in the multivariable MR analyses.

The concentration and content of medium HDL particles showed protective effects in univariable and multivariable MR analyses. In particular, adjusting for the traditional lipid risk factors did not attenuate the effect of traits related to medium HDL. The concentration of and total lipid in small HDL particles showed protective effects in multivariable MR analyses, though the effect sizes were smaller than those of the medium HDL traits. The mean diameter of HDL particles (HDL-D) had almost no effect on MI in the univariable MR analysis, but after adjusting for the traditional lipid risk factors, it showed a harmful effect.

Table 2 reports the estimated effects of M-HDL-P, S-HDL-P, HDL-D, and traditional lipid traits (TG, LDL-C, HDL-C, ApoB, ApoA1) in the multivariable MR analyses. To better understand the role of HDL subfractions and particle sizes, we also included in the table the results of the multivariate MR analyses for the traditional lipid risk factors only. Those baseline analyses suggested that HDL-C/ApoA1 had a weak, non-significant protective effect on CAD, which is consistent with prior studies (Holmes et al., 2015; Wang et al., 2020). Adding S-HDL-P to the MR analysis did not substantially alter the estimated effects of the traditional lipid traits. However, when M-HDL-P or HDL-D was included in the model, the estimated effects of M-HDL-P and HDL-D changed substantially. In particular, when M-HDL-P was included in the multivariable MR analyses, HDL-C/ApoA1 showed a harmful effect on CAD. When HDL-D was included, HDL-C/ApoA1 showed a protective effect.

We identified four genetic variants that are associated with S-HDL-P, M-HDL-P, or HDL-D, not associated with LDL-C or ApoB, and associated with CAD: rs838880 (SCARB1), rs737337 (DOCK6), rs2943641 (IRS1), and rs6065904 (PLTP) (Figure 3). These SNP-cis gene pairs are also supported by examining expression quantitative trait loci (eQTL) in the tissue-specific GTEx data (Appendix 4). The first three variants were not associated with S-HDL-P. However, they had uniformly positive associations with M-HDL-P, L-HDL-P, XL-HDL-P, HDL-D, ApoA1, and HDL-C, and a negative association with CAD. The last variant rs6065904 had positive associations with S-HDL-P and M-HDL-P, negative associations with L-HDL-P, XL-HDL-P, HDL-D, negative but smaller associations with ApoA1 and HDL-C, and a negative association with CAD.

Figure 3

Download asset Open asset

We also investigated the effects of lipoprotein subfractions and particle sizes on MI/CAD using multiple GWAS datasets, MR designs and statistical methods. The results are provided in Appendix 3 and are generally in agreement with the primary results reported above. The diagnostic plots for S-HDL-P and M-HDL-P did not suggest evidence of violations of the instrument strength independent of direct effect (InSIDE) assumption (Bowden et al., 2015) made by RAPS and GRAPPLE (Appendix 4).

Two published GWAS of the human lipidome [Kettunen2016, Davis2017] measured lipoprotein subfractions and particle sizes using NMR spectroscopy. We investigated the 82 lipid and lipoprotein traits measured in these studies that are related to very-low-density lipoprotein (VLDL), LDL, and HDL subfractions and particle sizes. All the subfraction traits are named using three components separated by hyphen: the first indicates the size (XS, S, M, L, XL, XXL); the second indicates the category according to the lipoprotein density (VLDL, LDL, IDL, HDL); the third indicates the measurement (C for total cholesterol, CE for cholesterol esters, FC for free cholesterol, L for total lipids, P for particle concentration, PL for phospholipids, TG for triglycerides). A full list of lipid and lipoprotein traits used in our study can be found in Appendix 1—table 1 below.

We estimated the genetic correlation between lipoprotein subfractions, particle sizes, and traditional lipid risk factors using the LD score regression (Li et al., 2016). Appendix 2—figure 1–3 show the estimated genetic correlation matrix between selected traits using different datasets. Below the figures, Appendix 2—table 1 shows the estimated genetic correlations of the lipoprotein subbfractions with the traditional lipid risk factors using the Davis GWAS. The results in Appendix 2—table 1 were then used to screen the traits as described in Materials and methods.

Appendix 2—figure 1

Download asset Open asset

Appendix 2—figure 2

Download asset Open asset

Appendix 2—figure 3

Download asset Open asset

We implemented several Mendelian randomization (MR) designs and statistical methods to estimate the causal effect of lipoprotein subfractions and particles sizes on coronary artery disease. In general, we adopted the three-sample summary data MR design described in Zhao et al., 2019b, Wang et al., 2020 and we swapped the roles of the GWAS datasets whenever permitted by the statistical methods. More specifically, the statistical methods we used for univariable MR (RAPS, IVW, weighted median) require that the GWAS datasets for obtaining instruments, SNP effects on the exposure, and SNP effects on the outcome must have no overlapping sample. The multivariable MR method we used (GRAPPLE) allows the exposure and outcome GWAS to be dependent and estimates the proportion of overlapping sample. However, GRAPPLE still requires that the selection GWAS uses an non-overlapping sample.

The MR designs we implemented in this study are summarized in Appendix 3—table 1. We considered two ways of instrument selection for univariable MR. In ‘traditional selection’, the traditional lipid traits were used to select the instruments for the corresponding subfraction traits. That is, HDL-C was used to select SNPs for HDL subfractions and particle size, LDL-C for IDL and LDL subfractions and particle size, and TG for VLDL subfractions and particle size. This tends to select more instruments because the GWAS for traditional lipid traits had a larger sample size. In ‘subfraction selection’, the instrumental SNPs were selected for each lipoprotein subfraction and particle size using the same or closest trait in the selection GWAS. For example, if the exposure under investigation is S-HDL-L but it is not measured in the Davis GWAS (if it is used for selection), S-HDL-P is used instead for instrument selection.

For multivariable MR, we considered two models with different sets of exposures: TG, LDL-C, HDL-C, and the subfraction/particle size under investigation; TG, ApoB, ApoA1, and the subfraction/particle size under investigation. SNPs were selected as potential instruments if they were associated (p-value $\le {10}^{-4}$) with at least one of the four exposures. LD clumping was then used to obtain independent instruments, as described in Materials and Methods.

We briefly comment on the statistical methods used in univariable MR. All the three methods we used—RAPS, IVW, weighted median—require that the exposure GWAS and outcome GWAS have non-overlapping samples. RAPS and weighted median can provide consistent estimate of the causal effect even when some of the genetic variants are not valid instruments, provided that the direct effects of the genetic variants are independent of the strength of their associations with the exposure. The last condition is called the Instrument Strength Independent of Direct Effect (InSIDE) assumption in the MR literature [bowden2015mendelian]. RAPS is also robust to idiosyncratically large direct effect (Bowden et al., 2015). Because IVW and weighted median can be severely biased by weak instruments (Zhao et al., 2020), we only used them with the set of SNPs that have genome-wide significant association (p-value $\le 5\times {10}^{-8}$) with the exposure. In comparison, RAPS does not suffer from weak instrument bias and we used it with all the SNPs obtained by LD clumping without any p-value threshold.

Below, Appendix 3—figure 1 shows the MR results for the 27 lipoprotein measurements selected in phenotypic screening. Estimates that are statistically significant at a false discovery rate of 0.05 are shown in Figure 2 of the main paper. Appendix 3—table 2 shows the estimated effect of all the lipoprotein subfractions and particle sizes on myocardial infarction or coronary artery disease in various MR designs. Full results of the multivariable MR analyses, including the estimated effects of the traditional lipid risk factors, can be found in Appendix 3—tables 5 and 6. The results of the univariable MR analyses using IVW and weighted median estimators can be found in Appendix 3—tables 3 and 4.

Appendix 3—figure 1

Download asset Open asset

In the tables below, Red indicates p-value is significant (at level 0.05) after Bonferroni correction for all the results in the corresponding table and blue indicates p-value ≤ 0.05.

Here we provide the list of modified Cochran's Q-statistics for the multivariable MR analyses (Appendix 3—tables 7 and 8).

As mentioned above, RAPS is more robust against invalid instruments than other statistical methods for univariable MR, but it still needs the InSIDE assumption to be approximately satisfied. Zhao et al., 2019b described two diagnostic plots RAPS that checks whether there is clear evidence that the InSIDE assumption is violated. Here, we report these plots for HDL-C and M-HDL-P in different studies (Appendix 4—figures 1 and 2). Notice that a lack of evidence to falsify the InSIDE assumption does not mean that it is true.

We can further assess the validity of the InSIDE assumption for M-HDL-P and S-HDL-P but examining the associations of their genetic instruments with the traditional lipid risk factors and other subfraction traits. We meta-analyzed the summary results in the two lipidome GWAS (Davis and Kettunen) and obtained SNPs that are associated with S-HDL-P and M-HDL-P (p-value $\le 5\times {10}^{-8}$; the results are LD-clumped). The next two Tables show some information about these genetic markers and their associations with other traits (Appendix 4—table 1 and 2).

Appendix 4—figures 3 and 4 shows how adjusting for LDL-C and TG changes the effects of the selected SNPs for S-HDL-P and M-HDL-P on CAD. The adjusted effect on CAD is obtained by original effect on CAD – 0.45 * effect on LDL-C – 0.25 * effect on TG. After the adjustment, the associations of the genetic variants with CAD generally became closer to the fitted lines that correspond to the estimated effects of S-HDL-P and M-HDL-P.

Appendix 4—figure 3

Download asset Open asset

Appendix 4—figure 4

Download asset Open asset

Here we provide evidence of variant-gene associations from Quantatitive Trait Locus (QTL) analyses in the GTEx project (Appendix 4—table 3).