- Research
- Open access
- Published:
Implicating genes, pleiotropy, and sexual dimorphism at blood lipid loci through multi-ancestry meta-analysis
Genome Biology volume 23, Article number: 268 (2022)
Abstract
Background
Genetic variants within nearly 1000 loci are known to contribute to modulation of blood lipid levels. However, the biological pathways underlying these associations are frequently unknown, limiting understanding of these findings and hindering downstream translational efforts such as drug target discovery.
Results
To expand our understanding of the underlying biological pathways and mechanisms controlling blood lipid levels, we leverage a large multi-ancestry meta-analysis (N = 1,654,960) of blood lipids to prioritize putative causal genes for 2286 lipid associations using six gene prediction approaches. Using phenome-wide association (PheWAS) scans, we identify relationships of genetically predicted lipid levels to other diseases and conditions. We confirm known pleiotropic associations with cardiovascular phenotypes and determine novel associations, notably with cholelithiasis risk. We perform sex-stratified GWAS meta-analysis of lipid levels and show that 3–5% of autosomal lipid-associated loci demonstrate sex-biased effects. Finally, we report 21 novel lipid loci identified on the X chromosome. Many of the sex-biased autosomal and X chromosome lipid loci show pleiotropic associations with sex hormones, emphasizing the role of hormone regulation in lipid metabolism.
Conclusions
Taken together, our findings provide insights into the biological mechanisms through which associated variants lead to altered lipid levels and potentially cardiovascular disease risk.
Background
Abnormal blood lipid levels are a major cause of cardiovascular disease [1], the leading cause of morbidity and mortality worldwide [2]. Conventional blood lipid measures, low-density lipoprotein cholesterol (LDL-C), total cholesterol (TC), triglyceride (TG), high-density lipoprotein cholesterol (HDL-C), and nonHDL-C (TC – HDL-C), are commonly used in clinical practice to identify individuals at high risk for cardiovascular events. Several treatments for reducing LDL-C, including statins, ezetimibe, and PCSK9 inhibitors [3], also reduce the risk of developing cardiovascular disease.
Genome-wide association studies (GWAS) for blood lipids have identified nearly 1000 associated genetic loci to date [4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23], including our recent multi-ancestry GWAS meta-analysis in 1.65 M individuals [24]. The latter focused on the gains from the multi-ancestry meta-analysis relative to the single-ancestry results, in terms of number of loci, fine-mapping, and polygenic score (PGS) transferability. However, a challenge in the field is that the underlying gene and biological pathways is often unknown for GWAS loci. Within lipid GWAS, prior fine-mapping studies combined with functional follow-up have successfully identified causal genes with high confidence for only a handful of associated GWAS loci, including SORT1 [25], TM6SF2 [12], and ANGPTL3 [26], among others. Highly sophisticated methods are emerging to prioritize causal genes in well-powered GWAS studies, such as the Data-driven Expression-Prioritized Integration for Complex Traits [27] (DEPICT) and the Polygenic Priority Score [28] (PoPS), that take into account genome-wide properties of associated loci and larger sets of associated loci are beneficial. These methods can be combined with algorithms that integrate expression data such as transcriptome-wide association studies (TWAS) and comprehensive experimental data sets such as mouse gene knockouts. Gene sets enriched for causal genes will enhance our ability to unravel the biological pathways underlying these associations and there is growing interest in using a combination of gene prioritization methods to provide compelling evidence for putative causal genes [29].
In parallel, the linkage of electronic health records with genetic data in large-scale population studies and patient biobanks allows for the systematic exploration of pleiotropy of lipid-associated alleles. While blood lipid levels have a well-documented causal effect on cardiovascular disease based on genetic association studies validated by randomized controlled trials [30,31,32], genetic pleiotropic associations might exist for other conditions. Unraveling such pleiotropy may yield new biological insights by revealing previously unrecognized connections between blood lipids and both cardiovascular and non-cardiovascular diseases. Phenome-wide association scans (PheWAS) adopt an agnostic approach to test for pleiotropic associations between genetic factors and a wide range of phenotypes [33]. Such knowledge may allow for the identification of lipid levels as novel diagnostic biomarkers, the repurposing of drugs, and the prevention of adverse drug events [34].
Finally, given empirical sex differences in blood lipid distributions, sex-specific genetic associations may yield novel biological insights. Pre-menopausal females have lower levels of LDL-C than same-age males, and HDL-C levels are higher among females of all ages compared to males [35]. Lipid levels also show a greater estimated heritability in females compared with males [36], especially for LDL-C and TC (> 1.3-fold difference). Sexual dimorphism in lipid levels may be partly explained by X chromosome variants. Evidence from human X-linked abnormalities (like Turner or Klinefelter syndromes) suggests an important role of this chromosome in lipid metabolism [37]. This is further corroborated by the lipid and atherosclerosis profiles in the Four Core Genotypes mouse model [38], which comprises XX and XY gonadal males and XX and XY gonadal females. GWAS studies have traditionally understudied the X chromosome due to technical and analytical difficulties. A recent high coverage whole X chromosome sequencing study [39] prioritized CHRDL1 as a candidate causal lipid gene, suggesting with larger sample sizes we may be able to discover additional variation on the X chromosome associated with lipid levels.
In this study, we first prioritize genes at GWAS lipid loci through multiple in silico gene prediction algorithms and experimental data sources using the latest Global Lipids Genetics Consortium multi-ancestry meta-analysis [24]. We then identify novel disease associations related to lipid levels through PheWAS in two large biobanks using PGSs. Finally, we perform sex-stratified meta-analysis to compare the associations between males and females to identify genetic loci with sex-specific associations and GWAS meta-analysis of the X chromosome, to better understand lipid level differences between the sexes. Together, our results highlight biological mechanisms through which lipid-associated variants lead to altered lipid levels.
Results
Identifying functional genes in lipid-associated loci
In a GWAS meta-analysis of blood lipid levels from 1.65 million individuals (Additional file 1: Table S1) at 91 million genotyped and imputed genetic variants, we observed a total of 2286 genome-wide significant index variants associated with lipid levels at 923 loci (± 500 kb regions). This corresponded to 416 index variants associated with LDL-C, 539 with HDL-C, 461 with TG, 487 with TC, and 383 with nonHDL-C. Uniquely, we observed 1750 variants associated with one or more lipid levels [24] (Additional file 2: Table S2).
We employed six approaches to identify candidate functional genes for all 2286 lipid associations. Our prioritization approaches include four locus-specific methods that are based on local information around the indexed variant: (1) the closest gene to the index variant, (2) genes with lipid-associated protein-altering variants, (3) colocalized expression quantitative trait loci (eQTL) genes, and (4) nearby genes prioritized by transcript-wide association studies (TWAS). We also used two genome-wide methods that leveraged genome-wide similarities of features: (1) gene-level Polygenic Priority Score (PoPS) [28] and (2) Data-driven Expression-Prioritized Integration for Complex Traits (DEPICT) [27]. We further combined the two genome-wide methods with the locus-specific methods to increase the confidence in prioritized genes: (1) PoPS intersects with any locus-specific methods (PoPS +), and (2) DEPICT intersects with any locus-specific methods (DEPICT +) (Fig. 1). Since the genome-wide gene prioritization approaches can prioritize different genes for different lipid types at the same locus, we report the gene prioritization results for all 2286 lipid-variant associations (Additional file 2: Table S2, Additional file 3: Figure S1).
We took the genes prioritized by PoPS + and performed text mining to determine whether previous biological evidence supported these genes as playing a role in lipid levels (Additional file 4: Table S3, S4). PoPS + leverages both locus-specific and genome-wide genetic signals to increase confidence level in prioritized genes [28]. In total, 882 out of 2286 lipid associations were assigned to one potential causal gene based on PoPS + . We identified a group of 466 unique genes among the 882 lipid associations. We determined that 31 out of the 466 PoPS + genes have over 1000 lipid-related publications, 91 PoPS + genes have 100–999 lipid-related publications, 321 PoPS + genes have 1–99 lipid-related publications, and 23 PoPS + genes had no lipid-related publications retrieved by the text mining algorithm. These 23 genes could indicate novel genes related to lipid levels for future work or be due to incorrect gene prioritization for a small fraction of index variants. (Additional file 5: Table S4). We then randomly selected 466 genes from 18,383 protein-coding genes using by the PoPS as the reference group to estimate the number of lipid-related publications we would expect to see by chance. A Mann–Whitney U test showed that there was a significant difference (W = 52,353, p-value < 2.2 x 10−16) between the set of genes identified by PoPS + compared to the reference set of 466 genes (Additional file 6: Figure S2). The median count of lipid-related publications was 19 for the PoPS + gene set compared with 2 lipid-related publications for genes in the reference set.
We performed a comprehensive lookup of all PoPS + prioritized lipid genes in the Therapeutic Target Database 2022 [24] and found 2092 drugs targeting at least one of our 102 PoPS + prioritized lipid genes observed in the database (Additional file 7: Table S5). Among those 102 PoPS + genes, we identify known drug target genes including PCSK9 druggable as subtilisin/kexin type 9 inhibitor, HMGCR druggable as HMG-CoA reductase inhibitor, PDE3A druggable as phosphodiesterase 3A inhibitor (CILOSTAZOL), and NR1H4 as a bile acid receptor FXR agonist (URSODIOL). We also identify several other potential drug targets [24] such as LIPG (lipase G) and NR1H3 (nuclear receptor subfamily 1 group H member 3), with relevant lipid biology. LIPG has phospholipase and triglyceride lipase activities and is a primary determinant of plasma HDL levels. NR1H3 has an important role in the regulation of cholesterol homeostasis, regulating cholesterol uptake through MYLIP-dependent ubiquitination of LDLR, VLDLR, and LRP8 that could be targeted as an LXR-alpha modulator.
Effects of protein-altering lipid alleles with protective effects on CAD, T2D, and NAFLD
Coronary artery disease (CAD), type 2 diabetes (T2D), and non-alcoholic fatty liver disease (NAFLD) are typically characterized by dyslipidemias. We examined protein-altering alleles with favorable lipid profiles for their associations with CAD, T2D, and NAFLD to identify potential cardiovascular drug targets without off-target liver or diabetes effects. Of the 2286 lipid associations, we observed 166 coding index variants. Eighteen coding variants with a protective lipid effect also had a protective effect for CAD or T2D (p-value < 0.001) and the lipid results colocalized with the CAD or T2D results, as appropriate, with a posterior probability of a shared causal variant > 0.8 (Table 1 and Additional file 8: Table S6). Six of these twenty variants had protective effects for both CAD and T2D, while nine were protective for CAD and three were protective for T2D (Table 1). Additionally, 269 noncoding alleles with a protective lipid effect also had a protective effect for CAD or T2D (p < 0.001; Additional file 8: Table S6).
Driver tissues for lipid levels
We applied DESE (Driver-tissue Estimation by Selective Expression) [40] to estimate the driver tissues of five lipid traits using both gene-level and transcript-level eQTL summary statistics from GTEx v8 tissues [41]. We identified liver as the top-ranked tissue for HDL-C (gene-level p-value = 4.5 x 10−18, transcript-level p-value = 3.0 x 10−26), TC (gene-level p-value = 1.1 x 10−25, transcript-level p-value = 1.4 x 10−33), and nonHDL-C (gene-level p-value = 2.0 x 10−19, transcript-level p-value = 3.9 x 10−29) based on both gene-level and transcript-level selective expression (Additional file 9: Figure S3, Additional file 10: Table S7). For LDL-C, we identified the spleen as the top-ranked tissue using GTEx gene-level data (p-value = 8.3 x 10−20), while liver was ranked second (p-value = 4.8 x 10−17). However, when using GTEx transcript-level data, liver was the top-ranked tissue (p-value = 4.3 x 10−29) and second was whole blood (p-value = 4.3 x 10−20). The top tissue for TG according to both GTEx gene-level and transcript-level expression data was whole blood (gene-level p-value = 6.4 x 10−20, transcript-level p-value = 1.4 x 10−21). Spleen and liver were second according to GTEx gene-level and transcript-level expression data, respectively. The results were consistent with previous knowledge that the liver is a major tissue for lipid metabolism. Transcript-level selective expression provided more statistically significant results for the estimated driver tissues compared to the gene-level selective expression, as reported in the original [40].
Polygenic scores for lipid phenotypes and phenome-wide association scans
We have previously reported that a polygenic score (PGS) for LDL-C was most informative when generated from the multi-ancestry GWAS and that the multi-ancestry PGS performed equally well in European-ancestry Americans, African-ancestry Americans, and continental Africans [24]. Using a similar approach, we generated PGS for the other four lipid traits (“Methods”).
We next performed a phenome-wide association scan (PheWAS) for the multi-ancestry lipid PGS (LDL-C PGS previously reported [24]) to identify pleiotropic effects of lipids with other traits in the European subsets of the UK Biobank and the Million Veteran Program (MVP) cohorts. We compared the effect sizes from the PheWAS analysis between the UK Biobank and MVP per lipid PGS and observed a moderate correlation between the two datasets (Additional file 11: Figure S4). The correlation of the PGS effects on all phenotypes between the UK Biobank and MVP ranges from 0.12 for the HDL-C PGS to 0.39 for the TC PGS (Additional file 11: Figure S4). In general, correlations were stronger for the ICD-10-based phecodes (r2 of 0.42–0.52) compared to the biomarkers (r2 of 0.06–0.23) (Additional file 11: Figure S4), which may reflect differences in study populations due to varied environmental effects, prevalence of chronic health conditions, and sex distribution. Among PheWAS results with p-value ≤ 0.05 in the UK Biobank, the correlation was even higher for ICD-10-based phecodes (r2 of 0.52–0.76) but remained the same for the biomarkers (r2 of 0.07–0.22).
We meta-analyzed the results from the two cohorts to increase the power of the PheWAS, by matching ICD10-mapped phecodes and biomarkers. In the combined the UK Biobank-MVP PheWAS results, we detected 58 phenotypes associated with the LDL-C PGS at phenome-wide significance level (p-value ≤ 6.5 × 10−5, corrected for 773 phenotypes), 165 with the HDL-C PGS, 59 with the TC PGS, 166 with the TG PGS, and 78 with the nonHDL-C PGS (Fig. 2, Additional file 12: Table S8, Additional file 13: Figure S5, Additional file 14: Figure S6, Additional file 15: Figure S7, Additional file 16: Figure S8). As expected, multiple cardiovascular phenotypes related to atherosclerosis, including the expected coronary artery disease as well as aortic aneurysm and essential hypertension, were phenome-wide significantly associated with all five lipid PGSs, indicating increased risk of these diseases for individuals with genetically predicted increased LDL-C, TG, TC, or nonHDL-C or genetically predicted decreased HDL-C. A recent wide-ranging Mendelian randomization analysis confirmed the causal effect of circulating lipids, not only for coronary artery disease, but other cardiovascular outcomes [42]. Additionally, all lipid PGSs were also significantly associated with decreased levels of direct bilirubin (Additional file 12: Table S8, Fig. 2, Additional file 13: Figure S5, Additional file 14: Figure S6, Additional file 15: Figure S7, Additional file 16: Figure S8), indicating genetically predicted lower LDL-C increased levels of bilirubin (Fig. 2). Correspondingly, lipid PGSs were associated with lower risk for cholelithiasis (gallstones) with the opposite direction for TG PGS, indicating that extreme lowering of LDL-C may impact rates of cholelithiasis (Additional file 12: Table S8, Fig. 2, Additional file 13: Figure S5, Additional file 14: Figure S6, Additional file 15: Figure S7, Additional file 16: Figure S8). To further clarify whether this association might be driven by the ABCG8 gene alone, we excluded from the LDL-PGS all variants within the locus and tested the association between LDL-PGS and cholelithiasis in the UK Biobank. There was no attenuation of the observed association (OR = 0.94 and p-value = 7.94 × 10−17 without the ABCG8 locus vs. OR = 0.93 and p-value = 1.96 × 10−21).
In the PheWAS analysis, we found that the TC and LDL-C PGS were significantly associated with increased levels of HbA1c (beta = 0.101 mmol/mol per SD PGS increase, p -value= 1.21 × 10−23 and beta = 0.095 mmol/mol per SD PGS increase, p-value = 4.37 × 10−21, respectively), while the HDL-C PGS was associated with decreased levels of HbA1c (beta = − 0.257 mmol/mol per SD PGS increase, p-value = 2.84 × 10−143) (Additional file 12: Table S8). Furthermore, genetically predicted increased LDL-C was significantly associated with decreased hemoglobin concentration (p-value = 1.92 × 10−45, similar significant associations for all other lipid PGSs with a reverse direction of effect for TG, Additional file 12: Table S8). As expected, genetically predicted increased LDL-C and TC were both associated with increased risk for Alzheimer’s disease [43] (OR = 1.33 per SD PGS increase, p-value = 1.74 × 10−44 and OR = 1.26 per SD PGS increase, p-value = 1.48 × 10−30, respectively; Additional file 12: Table S8). To further investigate how this association might be driven by the ApoE locus, we excluded all genetic variants overlapping this gene from the LDL-PGS and repeated the analysis in the UK Biobank. While the association between the LDL-PGS and the risk for Alzheimer’s disease was slightly attenuated after removing the ApoE locus (OR = 1.23 vs. 1.36 per SD PGS increase), the association remained significant (p-value = 2.51 × 10−21). Recent Mendelian randomization studies also provide evidence for the causal effect of lipids on risk for dementia [44] and Alzheimer’s disease [45]. The LDL-C and TC PGSs were also associated with increased aspartate aminotransferase levels (a liver enzyme), in accordance with other studies [46]. We also observed inverse associations between LDL-PGS (p-value = 1.43 × 10−14) and TC PGS (p-value = 8.34 × 10−14) with the risk of iron metabolism disorders (Additional file 12: Table S8).
To better understand the loci that drive the association between each of the lipid PGSs and cholelithiasis and cholecystitis, we interrogated the results from the single-variant PheWAS meta-analysis in the UK Biobank and MVP with all lipid multi-ancestry index variants (N = 1750 unique). We identified 22 genetic variants associated with cholelithiasis and/or cholecystitis at genome-wide significance. Genes prioritized for these index variants included genes already reported to be associated with gallstone disease [47] (CYP7A1, ABCG5/8), as well as additional genes (ABCB4, LRBA, HNF4A, NUCB1, GATA4), that may play also a role. Importantly, we found there was overlap (same index variant) between the previously published index variants for gallstone disease and our lipid index variants for these two loci (Additional file 17: Table S9).
Lipid loci show sex-specific effects
Sex-stratified analyses have the potential to identify loci missed by sex-combined analyses [48] as well as to detect loci exhibiting differential effects on lipids between sexes. First, we performed GWAS meta-analysis separately in each sex (Nmales = 749,391; Nfemales = 562,410), excluding loci discovered in the sex-combined analysis [24]. We identified twelve loci in females and four in males that reached genome-wide significance in the sex-stratified analysis (p-value < 5 × 10−8; Additional file 18: Table S10, Additional file 19: Table S11, Additional file 20: Table S12) but not in the sex-combined meta-analysis. As variants may show association to a single sex for reasons unrelated to biological sex differences, including differences in sample sizes between groups, we additionally tested for heterogeneity by sex for these variants in GLGC participating cohorts with close to equal number of males and females. Of the sixteen loci, eight showed significant sex-heterogeneity (p-value < 0.0031, Bonferroni-corrected threshold for sixteen tests). For example, the non-synonymous variant rs34372369 (EPHA1, p.Pro582Leu) is associated with nonHDL-C only in females (male p-value > 0.05) and shows significant sex-heterogeneity (p-value = 0.0012). This variant has been previously found to be linked with expression levels of the sex hormone-binding globulin gene (SHBG) more strongly in males than females [49], suggesting a possible reason for the difference in observed associations. We additionally sought to replicate the sex-heterogeneity results of these variants in 8 independent multi-ancestry cohorts (N = 311,639, 77% non-European ancestry, Additional file 21: Table S13). However, we did not detect significant differences in effect sizes between sexes for these variants after accounting for the number of tests (p-value > 0.0031, Additional file 22: Table S14), potentially due to the limited sample size or difference in ancestry makeup.
Second, we tested for a difference in the male- and female-specific effect sizes for each of the index variants identified from the sex-combined multi-ancestry meta-analysis. Of the 1750 unique index variants, 64 showed a significant difference in effect size by sex for one or more traits (Bonferroni correction for the number of index variants in each trait, Additional file 23: Table S15). These were evenly distributed across traits and more often had stronger effects in females than males (67%, Additional file 24: Figure S9). We tested for replication of the sex-specific differences in up to 311,120 participants from eight independent multi-ancestry cohorts not included in the original meta-analysis (Additional file 21: Table S13). Fifty-four of the 64 (84%) variants had effect size differences that were directionally consistent with the original meta-analysis (Additional file 25: Table S16). Of these, 10 had significantly different effect sizes (p-value < 7.8 × 10−4, Bonferroni correction for 64 variants) and 22 were nominally significant (p-value < 0.05). We attribute the low rate of replication to the small sample size and the differing proportions of ancestry groups within our replication samples, but we cannot dismiss the potential of false positives in the sex-specific discovery results.
We tested whether the observed sex differences could be caused by a higher frequency of cholesterol-lowering medications in males, potentially indicating an insufficient correction for pre-medication cholesterol levels. Among white British individuals in the UK Biobank, variants with significant sex differences had significantly higher effect size estimates on average after excluding individuals on medication (Additional file 26: Figure S10, Additional file 27: Table S17). However, of the 17 variants that exhibited a significant difference in effect size by sex in the UK Biobank alone, 15 remained significant after excluding individuals taking medications. Based on this observation, the observed differences did not appear to be driven solely, or even primarily, by differences in medication use between sexes. Furthermore, none of the identified sex-specific variants were associated with sex-participation bias [50] (Additional file 28: Table S18), indicating that differential study enrollment between sexes was unlikely to be the cause of the observed sex-specific lipid associations. We next investigated differences in environmental factors between sexes for these variants in the UK Biobank (Additional file 29: Table S19), including alcohol use [48], smoking status [48], body mass index (BMI) [51], and waist-hip ratio adjusted for BMI (WHRadjBMI) [51]. Twenty-two of the variants (34%) with differential effects on lipids by sex also exhibited a significant difference by sex for WHRadjBMI and one variant had a significant difference by sex for alcohol use (ADH1B p.His48Arg). The observed sex differences may therefore be partially attributed to pleiotropic associations with other traits.
Finally, we annotated each locus that showed significant sex differences with regulatory information to identify biological mechanisms that could underlie this difference. Of the 64 lipid variants with significant sex-stratified associations, 14 colocalized (posterior probability of H4 > 0.8) with expression of 20 genes in lipid-related tissues (liver, adipose, or skeletal muscle; Additional file 30: Table S20). Eight of these loci also show a sex-biased eQTL effect in at least one tissue in the direction concordant with the observed sex specificity of the GWAS effect (Additional file 30: Table S20). Among these ten is CETP, a gene with strong prior evidence for association with lipids, and UGT2B17 [20] (Additional file 31: Supplementary Note, Fig. 3). The lead variant of UGT2B17, rs4860987, shows a significantly stronger effect of LDL-C in males (Betamale = 0.042, SEmale = 0.002, Betafemale = 0.016, SEfemale = 0.003, p-valuedifference = 4.2 × 10−15) and colocalizes with a male-specific liver eQTL associated with increased expression of UGT2B17. Common variants at this locus are in moderate LD (R2 = 0.51) with a common copy number variation (CNV), which may mediate the causal effect (Additional file 31: Supplementary Note). UGT2B17 plays a role in the metabolism of androgens [52], including testosterone, which is consistent with the observed pleiotropic relationship of this locus with testosterone in males (Additional file 30: Table S20). We note that the index variant in UGT2B17, rs4860987, did not show significant sex-specific effects in the replication cohorts, but this could be due to varying frequencies for the index variant between ancestry groups and the moderate LD to the causal CNV in the region. We observed that the combined frequency of rs4860987 across the replication studies was much lower (8%) compared with our combined frequency in the discovery (24%) due to differing proportions of ancestry groups and, along with the lower number of individuals (N = 218,437), led to a much-reduced power to replicate this sex-specific effect.
Lipid-associated loci on the X chromosome
Lastly, we meta-analyzed association statistics for 3.1 million X chromosomal variants, including PAR regions, across 1,238,180 individuals from multiple ancestry groups. We identified 28 index variants significantly associated with lipid levels (Additional file 32: Table S21), of which 21 have not been previously reported [20, 39, 53] (15 for HDL-C, 4 for LDL-C, 4 for TC, 5 for TG and 4 for nonHDL-C, Table 2). Among these 28 loci, two have index variants with a minor allele frequency (MAF) < 1% and three index variants are missense mutations (in genes ARSL, TSPAN6, and G6PD), all of which are novel. We validated the identified X chromosomal associations in up to 255,475 individuals from seven multi-ancestry cohorts (Additional file 21: Table S13). Twenty index variants were at least nominally associated (p-value < 0.05), with five reaching genome-wide significance in the replication cohorts alone (p-value < 5 × 10−8, Additional file 32: Table S21).
We additionally considered potential sex differences for the X chromosome variants. A missense variant in RENBP with MAF = 2.5% reached genome-wide significance only in males but was not significant in the sex-combined meta-analysis or in the female-only analysis (p-value = 4.59 × 10−8, 0.003 and 0.2, respectively). We also observe three X chromosome loci with significant heterogeneity in effects between sexes; however, these were not significant in the replication cohorts alone, possibly due to the lower sample size (Bonferroni correction for the number of index variants in each trait, Additional file 32: Table S21).
Using a PheWAS approach in the UK Biobank, we found four of the novel loci to have pleiotropic associations with body composition traits (FAM9B [HDL-C], EDA2R [HDL, TG], TSPAN6 [LDL-C, TC], and DOCK11 [HDL-C]), four variants with coronary atherosclerosis and ischemic heart disease, three with immune-related biomarkers (SLC9A7 [HDL-C], CLCN5 [HDL-C], THOC2 [HDL-C]), and two with blood clotting-related biomarkers (KLF8 [TG], TEX11 [HDL-C]) (Additional file 32: Table S21). Interestingly, two of the three sex-biased X chromosome variants demonstrate the most significant association with testosterone of all lipid X chromosome variants tested in the PheWAS (rs505520: beta/SE = − 0.089/0.007 nmol/L per TG-increasing allele and rs5934507: beta/SE = 0.237/0.006 nmol/L per HDL-increasing allele).
Discussion
In this study, we identify and prioritize likely candidate genes at lipid-associated loci discovered through a variety of approaches including multi-ancestry meta-analysis of autosomes [24] (~ 91 million variants) and the X chromosome (~ 3 million variants), as well as sex-specific meta-analyses using sample sizes ranging from 1.35 to 1.65 million individuals. We previously reported a comparison of multi-ancestry vs single-ancestry lipid findings using autosomal chromosomes and identified improvements in fine-mapping of credible sets and PGS performance, with slight differences in the number of identified loci by ancestry group [24]. Here, we add X chromosome and sex-specific discovery results. We also focus on lipid biology by prioritizing implicated genes, identifying novel phenotypes and diseases associated with genetically predicted lipid levels, and predicting candidate drug target genes.
Our results from this effort translate our GWAS findings for three complimentary research areas, helping us further elucidate the biological mechanisms underlying the lipid-associated genetic variants. We first sought to identify methods for prioritization of functional genes at GWAS loci by performing six gene prioritization methods. Lipids are an excellent exemplar phenotype for gene prioritization algorithms because of a wealth of GWAS loci (~ 1000), Mendelian dyslipidemia genes (21), and mouse dyslipidemia phenotypes observed in gene knockouts (740). While the gene prioritization approaches are not independent of each other, integrating several prioritization predictors provides higher confidence when attempting to characterize causal genes. Others have also highlighted the importance of such frameworks in different diseases [29, 54, 55].
We identify 466 unique genes by combining evidence from a global approach (PoPS) with local gene prioritization approaches. The vast majority of these genes had many lipid-related publications, suggesting the accuracy of our combined prioritization approach. Twenty-three PoPS + identified genes had no lipid-related publications, indicating they could be truly novel or possibly were incorrectly prioritized. Functional validation of the larger pool of prioritized genes, which will require highly parallel experimental methods, will help to further optimize bioinformatics algorithms to prioritize genes and is beyond the scope of this manuscript.
Our prioritization approach also indicates several genes as potential drug targets including PDE3A and NR1H4. PDE3A encodes the phosphodiesterase 3A gene and is predicted to be druggable as phosphodiesterase 3A inhibitor (CILOSTAZOL). Cilostazol has antiplatelet, anti-proliferative, vasodilatory, and ischemic-reperfusion protective properties [56] and has been previously suggested for the primary or secondary prevention of CAD [22]. NR1H4 encodes a bile acid receptor and regulates the expression of genes involved in bile acid synthesis and transport. The target gene is predicted to be druggable as a bile acid receptor FXR agonist (URSODIOL). Ursodiol is used to treat primary biliary cirrhosis and cholelithiasis and could be a potential candidate for drug repurposing.
We also identify eighteen coding variants where the protective lipid allele is also protective for CAD or T2D. Among these, PCSK9 is a well-documented drug target, not only for lipids but also for cardiovascular events [57,58,59]. In comparison to published studies [60], others find a non-significant increased risk for T2D [61] and an arguably stronger protective effect for CAD [62], for PCSK9 variant carriers. Our observation is consistent with the lack of excess T2D risk observed in PCSK9 inhibitor clinical trials [57,58,59, 63] and with strong protective effects for coronary heart disease [64]. Furthermore, these variants are potential therapeutic targets for protective lipid profiles and lowering risk of disease.
Our second goal was to identify diseases that may benefit from lipid-lowering as well as diseases or traits that may become problematic due to very low lipids. To accomplish this, we examined the association of genetically predicted lipid traits (using PGS) with 773 phenotypes in 478,556 individuals. We observed that genetically predicted increased LDL-C, TC, and HDL-C levels, or decreased TG levels, decrease the risk of cholelithiasis. Prior epidemiological studies have not consistently reported an association between lipid levels and risk of gallstones, with some studies showing that increased levels of LDL-C, TC, and TG and decreased levels of HDL-C predispose to the risk for cholelithiasis [65, 66], but others showing no association [67, 68]. Our results are corroborated by a recent Mendelian randomization meta-analysis study in the FinGen and UK Biobank cohorts [69]. The prioritized genes for the individual index lipid variants significantly associated with cholelithiasis in the PheWAS analysis include ABCG8, a hepatic cholesterol transporter, responsible for the efflux of cholesterol from the enterocytes to the lumen and from the hepatocytes into bile [70]. The lipid-decreasing allele of index variant in ABCG8, rs4245791, has been previously associated with high risk for gallstone disease [47] and high intestinal cholesterol absorption [71], possibly mediated by an increased expression of ABCG8 [72]. Furthermore, even after excluding ABCG5/8 variants from the LDL-PGS, the association with the risk of cholelithiasis was not attenuated. These PGS-PheWAS results suggest the existence of many other cholesterol transporters like ABCG8 that modify blood cholesterol levels perhaps in large part by facilitating an increased secretion of cholesterol into the biliary system, which in turn increases the risk of the formation of gallstones through the supersaturation of bile. We also observed that HbA1c levels were elevated among subjects with genetically predicted increased LDL-C and TC and with genetically predicted decreased HDL-C. Previous epidemiological studies have established associations between dyslipidemia (increased LDL-C, TC, TG, and decreased HDL-C levels) and increased HbA1c levels among subjects with T2D, as well as insulin-resistant subjects without diabetes [73, 74]. Our observations support a strong genetic basis to these associations and are in accordance with previous studies showing shared pathways between lipid biology, T2D, and HbA1c [75], as well as pleiotropic effects of blood red cell variants with lipid levels [76]. Mendelian randomization studies have shown that hemoglobin and LDL show bidirectional inverse relationships and hemoglobin effects on LDL are also mediated through Hb1Ac, implying that genetic variation influencing erythrocytic factors could also determine lipid levels and the opposite [77]. While most of our significant PheWAS findings could be confirmed via Mendelian randomization studies, we cannot exclude the possibility of spurious associations due to pleiotropy.
Lastly, we sought to expand the coverage of the genome and performed the most comprehensive GWAS of lipid levels to date by including assessment of 3 million variants on the X chromosome as well as explicitly testing for sex-specific effects across 23 chromosomes in 1.35 million individuals of diverse ancestries. We report 21 novel X chromosome loci, including an LDL-lowering locus involving a missense variant in G6PD, known to cause G6PD deficiency (p.V68M) [78]. The proposed mechanism is via the inhibition of the NADPH-dependent hydroxymethylglutaryl-CoA (HMG-CoA) reductase, resulting in decreased cholesterol biosynthesis, even though the protective effect of the G6PD deficiency on cardiovascular risk is debatable [79].
We also observed that approximately 3–5% of the genome-wide lipid index variants exhibited differential effects between sexes, which did not seem to be due to differential prevalence in the use of lipid medications or study selection bias. These findings may have important implications in the interpretation of lipid biology, the identification of novel drug targets, and possibly for more accurate prediction of blood cholesterol-related conditions. For example, the UGT2B17 locus, one of the ten sex-biased loci with corresponding sex-biased eQTL effect, is known to be implicated in androgen and drug metabolism [52]. A common CNV in the region, partially tagged by the lipid index variant, is associated with significant variations in expression levels between ethnic groups [80], which would explain lack of replication in the set of independent studies, and the deletion has been linked to testosterone-related decreased BMI levels [81], as well as decreased risk for osteoporosis in men [82].
Several of the reported sex-biased and X chromosome loci showed significant pleiotropic effects with sex hormone levels, including testosterone and SHBG, highlighting the role of hormone regulation in lipid metabolism [83]. In particular, we observe an overall inverse effect between the X chromosome lipid index variants and the sex hormone levels. Observational studies have long suggested a potential influence of the sex hormones on the risk for cardiovascular risk [84] but this hypothesis has not been consistently supported by recent Mendelian randomization studies, raising the issues of reverse causality [85, 86].
Conclusions
In conclusion, we leverage the power of a large multi-ancestry GWAS study to further our understanding of lipid metabolism and the impact on chronic diseases. We identify novel lipid loci on the X chromosome and autosomal loci with evident sex-biased lipid effects. We compare a range of gene prioritizing methods to identify causal genes, an approach applicable to studying other complex traits. We additionally further our understanding of lipid metabolism through a phenome-wide study that implicates a relationship between genetically predicted low cholesterol with risk of cholelithiasis.
Methods
Meta-analysis
Summary statistics for sex-combined autosomal analyses were previously published [24]. Following the same procedure, we carried out meta-analyses stratified by sex for 5 lipid traits (HDL-C, LDL-C, TG, nonHDL-C, and TC) for both the autosomes and chromosome X. The sample size for chromosome X (Total N = 1,238,180; males = 749,391; females = 562,410) was lower than available for autosomes as not all participating biobanks submitted results for chromosome X. Quality control of summary statistics from contributing cohorts was performed using EasyQC [87]. Prior to meta-analysis, we removed variants with low imputation info scores (r2 < 0.3), those with minor allele count < 3, and those with Hardy–Weinberg equilibrium p-value < 1 × 10−8. Variants on the X chromosome were filtered using the female imputation info scores and Hardy–Weinberg equilibrium p-values. Summary statistics were corrected by the genomic-control factor calculated from the median p-value of variants with minor allele frequency > 0.5%. For cohorts that contributed summary statistics imputed both on the Haplotype Reference Consortium (HRC) and 1000 Genomes Population v3 (1KGP3) panels, we generated a single file containing all possible variants, favoring those imputed from the HRC imputation panel due to generally higher imputation quality of these variants. Multi-ancestry meta-analysis was performed with MR-MEGA [88] with 5 principal components and using the inverse-variance weighted method in METAL to estimate effect sizes [89]. Independent loci were defined with physical distance > 500 kb or genetic distance > 0.25 cM, whichever one would result in a larger window, followed by a conditional analysis using rareGWAMA [90] as previously described [24], to identify index variants that were shadows of nearby, more-significant associations. Conditional analysis for chromosome X used a female-only UK Biobank LD reference (N = 21,510). In line with the analysis in the autosomes, a locus was identified as dependent if the effect size after conditioning on the most significant variant in the area was more than 1.43 times smaller than the original (95th percentile of the effect size ratios for chromosome X).
Differences in effect size between males and females were tested within each cohort using [91]:
and were then meta-analyzed across studies using METAL, to account for cohort-specific ascertainment (e.g., enrichment of cases for type 2 diabetes), or demographics, such as age.
Replication
We collected summary statistics from 8 cohorts across 6 ancestry groups, including African or African American, East Asian, European, Hispanic, Middle Eastern, and South Asian. Each cohort provided sex-stratified and X chromosome association results for the tested traits, as available. The difference in effect sizes between males and females was calculated within each cohort as described above and then meta-analyzed across studies using METAL. X chromosome association results were meta-analyzed using METAL with weighting by sample size.
Gene prioritization methods
Closest gene
We annotated the closest gene to the lipid multi-ancestry index variants [24] by identifying the closest gene transcript on either side (500 kb) of the index variant [92].
Colocalization with GTEx eQTLs
For each of the five lipid phenotypes, we first lifted over GWAS summary statistics from the multi-ancestry meta-analysis [24] to GRCh38 using the UCSC liftOver tool. Then, we defined a set of approximately independent windows across the genome within which colocalization with eQTLs was run. To define these, we first identified all genome-wide significant variants (p-value < 5e − 08) from the meta-analysis for each lipid trait and sorted them by significance, from most significant to least. Starting with the most significant variant, we aimed to define a window defining independent genetic signals; we define a variant’s window as a region within the greater of 500 kb or 0.25 cM on either side of this “sentinel variant.” Genetic distances were defined using reference maps from HapMap 3. All other genome-wide significant variants within this window were discarded from the list of sentinel variants, and similar windows were defined for the remaining genome-wide significant variants.
We ran an eQTL colocalization using GTEx v8 eQTL summary statistics within each of our defined windows for all lipid traits. For each of the 49 GTEx tissues, we first identified all genes within 1 Mb of the sentinel variant, and then restricted analysis to those genes with eQTLs (“eGenes”) in that tissue (FDR < 0.05). We used the R package “coloc” (run on R version 3.4.3, coloc version 3.2.1) [93] with default parameters to run colocalization between the GWAS signal and the eQTL signal for each of these cis-eGenes, using as input those variants in the defined window, i.e., all variants present in both datasets. A colocalization posterior probability of (PP3 + PP4) > 0.8 was used to identify loci with enough colocalization power, and PP4/PP3 > 0.9 was used to define those loci that show significant colocalization [94].
Transcriptome-wide association studies (TWAS)
For our transcriptome-wide association analysis (TWAS), we integrated the results of our GWAS with eQTL summary statistics from GTEx v8. The S-PrediXcan software [95] allows us to integrate these two datasets using only summary statistics from GWAS without needing individual-level genotype data. We used the multi-ancestry lipid GWAS summary statistics [24] and harmonized them with the GTEx summary statistics. Then we performed the TWAS using the eQTL models estimated on GTEx v8 expression data. For each of the 49 GTEx tissues, we identified “significant genes” those genes with p-values more significant than an FDR threshold of 0.05.
Genes with coding variants
We determine the coding variants within 99% credible sets and the coding variants in LD > 0.8 with variants in the 99% credible sets with the credible sets as defined here [24]. We define regions for construction of the credible sets as ± 500 kb around each index variant. We used Bayes factors (BFs) for each variant from the MR-MEGA output and generated the credible sets within each region by ranking all variants by BF and calculating the number of variants required to reach a cumulative probability of at least 99%. Additionally, we used previously established gene-based associations [96] to determine whether rare coding variation in a gene were associated with blood lipid levels (p < 0.001). We labeled a gene as having coding variants if any of these criteria were met.
DEPICT
We used Data-driven Expression-Prioritized Integration for Complex Traits (DEPICT, v1 beta version rel194 for 1 KG imputed GWAS) to prioritize genes at our index variants, on the assumption that truly associated genes share functional annotations [27]. Index variants [24] with p-value < 5 x 10−8 were retained as input. We implemented the DEPICT analysis with the default settings of 500 permutations for bias adjustment and 20 replications for FDR estimation. DEPICT prioritizes genes by calculating the similarity of a given gene to genes from other associated loci across 14,461 reconstituted gene sets and estimates the nominal gene prioritization p-value and the estimated false discovery rate of each gene. The prioritized genes at FDR < 0.05 were considered significant.
PoPS
We used the PoPS method to prioritize genes in the previously reported [24] multi-ancestry index variants for all lipid traits. The PoPS method [28] is a new gene prioritization method that identifies the causal genes by integrating GWAS summary statistics with gene expression, biological pathway, and predicted protein–protein interaction data. First, as part of the PoPS analysis, we used MAGMA to compute gene association statistics (z-scores) and gene–gene correlations from GWAS summary statistics and LD information from a multi-ancestry reference panel (1000 Genomes). Next, PoPS performs marginal feature selection by using MAGMA to perform enrichment analysis for each gene feature separately. The model is fitted by generalized least squares (GLS), and MAGMA results are used to perform marginal feature selection, retaining only features that pass a nominal significance threshold (p < 0.05). Then PoPS computes a joint enrichment of all selected features simultaneously in a leave one chromosome out (LOCO) framework. The gene features employed by PoPS are listed here: https://github.com/FinucaneLab/gene_features. Finally, PoPS computes polygenic priority scores for each gene by fitting a joint model for the enrichment of all selected features. The PoPS score for a gene is independent of the GWAS data on the chromosome where the gene is located. The PoPS analysis returned scores for a total of 18,383 genes per lipid trait. We only kept the top 20% genes among all 18,383 genes. We then annotated our index variants with the nearest ENSEMBL genes in a 500-kb window (either side) and selected the highest PoPS score gene in the locus as the prioritized one.
We performed the PoPS analysis on our lipid-specific multi-ancestry meta-analysis results, using all populations from 1000G as the reference for the LD information in MAGMA. As a sensitivity step, we also repeated the same analysis using only the European population from 1000G as the reference. We observed high concordance in the top two PoPS prioritized genes from both reference panels. In detail, the same 2119 genes (89%) were prioritized as the top genes from both panels, a further 203 genes were prioritized as a top gene with one panel and as the second top with the other and only 7 genes were completely mismatched between the two reference panels.
Monogenic genes
We annotated genes known to cause Mendelian lipid disorders based on proximity with identified GWAS loci [97, 98]. GWAS index variants within ± 500 kb of the transcription start and end positions from the USCS genome browser annotations were annotated as nearby known monogenic dyslipidemia genes.
Mouse knockout lipid phenotype silver set genes
Human gene symbols (9557 unique genes) were mapped to gene identifiers (HGNC) and their corresponding mouse ortholog genes were obtained using Ensembl (www.ensembl.org). Phenotype data for single-gene knockout mouse models were obtained from the International Mouse Phenotyping Consortium (IMPC) (www.mousephenotype.org) latest data release 12.0 (www.mousephenotype.org/data/release). The knockout mouse models were primarily produced by IMPC but also include some models that have been reported from the relevant literature and were curated by Mouse Genome Informatics (MGI) data release 6.16 (www.informatics.jax.org). For each mouse model, reported phenotypes were grouped using the mammalian phenotype ontology hierarchy into broad categories relevant to lipids: growth and body weight (MP:0001259), lipid homeostasis (MP:0002118), cholesterol homeostasis (MP:0005278), and lipid metabolism (MP:0013245). This resulted in mapping of human genes to significant phenotypes in animals.
For each of the multi-ancestry lipid index variant [24], we mapped the closest gene to the knockout mouse phenotypes and curated the set to only include mouse phenotypes strictly relating to lipid metabolism. That resulted in our silver set of 740 genes with mouse lipid phenotypes (Additional file 33: Table S22).
Overlap between methods
We standardized the gene names across different methods using the R/geneSynonym package, a wrapper to gene synonym information in ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene_info.gz. We also quantified how often the same gene was prioritized by multiple methods for each index variant and determined scores that ranged from 1 to 6 (S1-S6), based on the number of methods that prioritized the gene.
We integrated multiple gene prioritization methods to identify likely causal genes in the latest global lipid genetics consortium GWAS results. In total, we have implemented the 6 individual gene prioritization methods above that utilize the GWAS summary statistics from meta-analysis. Our gene prioritization methods can be placed into two broad categories, the locus-specific methods and the genome-wide methods. The locus-specific methods leverage local GWAS data by connecting the causal variants to the causal gene(s) using genomic distance, eQTLs, or protein-coding variants.
More specifically, there are four locus-specific methods that have been implemented including: (1) The closest protein-coding gene around the index variants based on the genomic distance, (2) eQTL colocalization using r COLOC package, (3) TWAS using S-PrediXcan, (4) coding variants which have been identified in 99% credible sets OR in LD > 0.8 with coding variants OR from gene-based tests (p < 0.001) of rare coding variants. For the eQTL and TWAS, we first used all the 49 GTEx tissues and then restricted to only 5 lipid-specific tissues: liver, adipose subcutaneous, adipose visceral, whole blood, and small intestine. In addition, two genome-wide methods were employed: (1) DEPICT (FDR < 0.05), (2) PoPS (Top 1 gene). It is reasonable to combine similarity-based methods with locus-based methods since they use two different sources of information.
To determine the relative performance of each prioritization method and their combined scores for lipid loci, we used 21 genes known to cause Mendelian dyslipidemias as a gold standard set (ABCA1, ABCG5, ANGPTL3, APOA5, APOB, APOE, CETP, CYP27A1, GPD1, GPIHBP1, LCAT, LDLR, LDLRAP1, LIPA, LIPC, LMF1, LPL, MTTP, PCSK9, SAR1B, SCARB1), and 740 mouse knockout genes causing lipid phenotypes as a silver standard set (Additional file 33: Table S22). We examined two metrics for each gene prioritization approach: (1) the proportion of prioritized genes in the gold/silver standard set, and (2) the proportion of correctly identified genes among all prioritized genes (Additional file 3: Figure S1). Note that out of the 2286 lipid associations, 97 fell within 500 kb of a Mendelian gene and 1280 within 500 kb of a mouse knockout gene with a lipid phenotype. We observed that the TWAS results yielded a high number of prioritized genes, but lead to a low proportion correctly identified. The TWAS approach had a much smaller proportion of genes correctly prioritized among all the prioritized genes, given it prioritized a total of 3511 genes, which was 3.5-fold greater than the other methods (~ 1000 genes). Notably, PoPS provided a similar proportion of correctly identified genes (78%) as of TWAS, while retained a high proportion of prioritized genes in the gold standard set (67%). Compared with PoPS, PoPS + (PoPS plus one of the local methods) slightly sacrificed the proportion of correctly identified genes from 78 to 71%, but improved the proportion of prioritized genes in the gold standard set from 67% to 79%. Overall, PoPS/PoPS + outperform other gene prioritization methods on both metrics for our gold (Additional file 3: Figure S1A) and silver (Additional file 3: Figure S1B) standard gene sets. We also assessed lipid-relevant tissue (liver, subcutaneous and visceral adipose, whole blood, and small intestine) expression QTLs (lipid eQTLs) and transcriptome-wide association (lipid TWAS) and found that the expression results from all tissues performed slightly better at recovering the reference gene sets compared with limiting to the lipid-relevant tissues (Additional file 3: Figure S1).
Text mining analysis
We retrieved the whole MEDLINE/PubMed titles and abstracts as of March 06, 2022, from National Library of Medicine (https://ftp.ncbi.nlm.nih.gov/pubmed/baseline/; https://ftp.ncbi.nlm.nih.gov/pubmed/updatefiles/). We then examined whether a list of genes prioritized by PoPS + and any one of the lipid-related keywords (lipid, lipids, triglyceride, triglycerides, fatty acid, cholesterol, dyslipidemias, hyperlipidemia, hypercholesteremia, diabetes, type 2 diabetes, type II diabetes, heart, cardiovascular, artery, coronary, coronary artery, coronary heart, atherosclerosis, peripheral vascular, PAD, stroke) occurred in the same abstract. We counted how many lipid-related publications that have a specific gene co-occurred with at least one lipid-related keyword. The same text mining approach was also implemented to a set of randomly selected genes from the 18,383 protein-coding genes used by the PoPS. We estimated the number of lipid-related publications we would expect to see by chance. A Mann–Whitney U test was performed to show whether there was a significant difference between the number of lipid-related publications of the PoPS + gene set and reference gene set.
Drug target mining analysis
To gain therapeutic insights from our gene prioritization results, we performed a lookup in Therapeutic Target Database (TTD) 2022 [99] (http://db.idrblab.net/ttd/). Specifically, we cross-referenced 466 unique lipid-associated genes prioritized by PoPS + (Additional file 2: Table S2) with 1563 genes corresponding to at least one drug (either under development or approved) with known clinical indication in TTD 2022. As a quality control for this lookup, we excluded all TTD entries related to drugs that were discontinued, terminated, or withdrawn from the market. The full lookup results are available in Additional file 8: Table S6.
Driver tissues for lipid levels
We performed phenotype-tissue association analysis using DESE (driver-tissue estimation by selective expression) [40]. DESE estimates the causal tissues by selective expression of phenotype-associated genes in GWAS. We used the GWAS summary statistics from the five lipid traits and the GTEx v8 normalized gene-level and transcript-level expression datasets as input. SNPs inside a gene and its ± 5 kb adjacent regions were first mapped to the gene, and then DESE ran iteratively to produce a list of driver tissues and the corresponding p-values of the associations. We used a Bonferroni-corrected significance threshold of 0.05/54 = 9.3 x 10−4.
PheWAS analysis
Construction of lipid PGSs
We had previously developed a multi-ancestry PGS for LDL-C that was demonstrated to perform well across multiple ancestry groups [24]. In a similar manner, we also generated PGS for HDL-C, nonHDL-C, TC, and triglycerides. First, multi-ancestry meta-analysis results were generated with METAL [89] after excluding individuals from the Michigan Genomics Initiative and the UK Biobank. The set of variants used to construct the PGS was limited to those that were well-imputed (R2 > 0.3) in MGI, UK Biobank, and MVP. Risk scores based on PRS-CS [100] or pruning and thresholding with Plink [101] across several r2 (0.1, 0.2), distance (250 kb, 500 kb), and p-value thresholds (5 × 10−10, 5 × 10−9, 5 × 10−8, 5 × 10−7, 5 × 10−6, 5 × 10−5, 5 × 10−4, 5 × 10−3, 0.05) were developed. For each trait, the single best score was selected based on the adjusted r2 calculated in the UK Biobank of the linear model for the lipid trait with the risk score and age, sex, batch, and PC1-4 as covariates. This corresponded to PRS-CS for HDL-C and nonHDL-C and pruning and thresholding for LDL-C (r2 = 0.1, p-value = 5 × 10−4, 500 kb), TG (r2 = 0.1, p-value = 5 × 10−3, 500 kb), and TC (r2 = 0.1, p-value = 5 × 10−4, 500 kb). The variance explained by the risk score among the UK Biobank participants was similar across traits (adjusted r2 of the full model-adjusted r2 of covariates: HDL-C = 0.13; LDL-C = 0.15; nonHDL-C = 0.14; TC = 0.14; TG = 0.10) and validated the ability of the risk score to predict genetically increased lipid levels.
PheWAS of lipid PGSs and index lipid variants in the UK Biobank and MVP
We used the European ancestry subset of individuals from the UK Biobank (408,886 samples) and the European samples from MVP (69,670 samples) to perform the PheWAS analysis.
We constructed a weighted PGS for each of the lipid traits, based on the corresponding genome-wide significant multi-ancestry index variants. We used the PheWAS package in R [102] to map ICD-10 codes from hospital records into clinically relevant phenotypes (phecodes) and to implement these association analyses, while adjusting for sex, age, 10 genetic principal components, and genotyping array (for the UK Biobank only) in each cohort. For the lipid-PGS PheWAS, each PGS was inverse normalized prior to analysis and lipid levels were corrected for statin use. The MVP samples used for the PheWAS analysis were not included in the GWAS meta-analysis [24].
Similarly, we extracted all multi-ancestry autosomal index variants for all lipid traits from the same European ancestry subset of the UK Biobank and MVP and performed a single-variant PheWAS association analysis per cohort. Additionally, we performed a single-variant PheWAS association analysis in the UK Biobank only with the sex-stratified and X chromosome index variants from the multi-ancestry analysis.
Meta-analysis of MVP and the UK Biobank PheWAS results
We combined, via meta-analysis, PheWAS lipid-specific PGS results for all intersecting phecodes and biomarkers between the UK Biobank and MVP (Europeans only) per lipid trait. We used ICD10-based phecodes and manually matched biomarkers to identify intersecting phenotypes between the two datasets. We restricted our meta-analysis to phenotypes that had at least 100 samples (total number for continuous traits or number of cases for binary traits) in each cohort. After the meta-analysis, we excluded phenotypes that had less than 500 combined samples (total number for continuous traits or number of cases for binary traits), to avoid reporting spurious results [103]. That resulted in a total of 773 phenotypes (739 phecodes and 34 biomarkers/measurements). We used both fixed and random effects model for the meta-analysis. We assessed heterogeneity using the p-value for Cochran’s q and set the level for significant heterogeneity at a Bonferroni threshold (p-value ≤ 6.5 × 10−5, to account for multiple testing of 773 phenotypes). We report the results from the fixed-effects model for the phenotypes with non-significant heterogeneity and the results from the random effects model for all others. Similarly, we meta-analyzed all index-variant PheWAS results between the UK Biobank and MVP and obtained results for 811 phenotypes and 1750 lipid multi-ancestry index variants, after excluding instances with a combined sample size < 500.
Lipid index variants with CAD, T2D, and NAFLD datasets
The GWAS meta-analysis results of CAD and T2D were acquired from MVP [62] and DIAGRAM Consortium [61], respectively. For variant rs1229984, the CAD result is from CARDIoGRAMPlusC4D meta-analysis [104], as it was not present in the MVP results. The NAFLD GWAS and meta-analysis was performed in the UK Biobank and Michigan Genomics Initiative (MGI). We determined the association of the lipid index variants with CAD, T2D, and NAFLD and aligned the alleles across all the traits to the LDL-lowering allele. We then highlighted the protective lipid coding alleles associated with CAD.
GWAS and meta-analysis of NAFLD in the UK Biobank and Michigan Genomics Initiative (MGI)
Individuals with NAFLD were identified using ICD-9 571.8 and ICD-10 K76.0. Individuals with hepatitis, liver cirrhosis, liver abscess, ascites, a liver transplant, hepatomegaly, jaundice, or with abnormal result of serum enzyme levels or a function study of the liver were excluded (exclusion phecodes 70.2, 70.3, 571.51, 571.6, 571.8, 571.81, 572, 573, 573.2, 573.3, 573.5, 573.7, 573.9) [105]. Analysis was performed using SAIGE v43.3 [106]. Analysis in the UK Biobank included white British individuals with batch, sex, birth year, and the first 4 genetic principal components as covariates. A total of 1122 cases and 399,900 controls were included in the analysis. Analysis in MGI included only European-ancestry participants with array version, sex, birth year, and the first 4 genetic principal components as covariates. A total of 2901 cases and 49,098 controls were analyzed. Meta-analysis was performed using METAL with weighting based on the effective sample size calculated as 4/((1/Ncases) + (1/Ncontrols)).
CAD/T2D colocalization analysis with lipid traits
We used R package coloc v3.2.1 [93] to perform summary statistics-based colocalization via a Bayesian approach and test whether the 5 lipid traits share common genetic causal variants with CAD or T2D. We first defined a window of ± 100 kb around each index variant [24]. Then for each window of the 10 pairs of traits, we ran colocalization with default parameters using those SNPs present in both datasets. A colocalization posterior probability of PP4 > 0.8 was used to define those loci that show significant colocalization.
Availability of data and materials
The GWAS meta-analysis results (including both ancestry-specific and trans-ancestry analyses) and risk score weights are available at: http://csg.sph.umich.edu/willer/public/glgc-lipids2021 [107]. A web browser displaying the gene prioritization and PheWAS results is available at https://hugeamp.org:8000/research.html?pageid=GLGC_149 [108]. The optimized trans-ancestry polygenic score weights are deposited within the PGS Catalog (https://www.pgscatalog.org/publication/PGP000230/ [109] and https://www.pgscatalog.org/publication/PGP000366/ [110]. Scripts used for analysis and summary of results are available under the MIT license on this GitHub repository: https://github.com/Global-Lipids-Genetics [111]. The version of source code used in the manuscript is deposited in Zenodo: https://doi.org/10.5281/zenodo.7130299 [112].
References
Castelli WP, Anderson K, Wilson PW, Levy D. Lipids and risk of coronary heart disease. The Framingham Study. Ann Epidemiol. 1992;2:23–8.
GBD. Diseases and Injuries Collaborators: Global burden of 369 diseases and injuries in 204 countries and territories, 1990–2019: a systematic analysis for the Global Burden of Disease Study 2019. Lancet. 2019;2020(396):1204–22.
Grundy SM, Stone NJ, Bailey AL, Beam C, Birtcher KK, Blumenthal RS, Braun LT, de Ferranti S, Faiella-Tommasino J, Forman DE, et al. 2018 AHA/ACC/AACVPR/AAPA/ABC/ACPM/ADA/AGS/APhA/ASPC/NLA/PCNA Guideline on the Management of Blood Cholesterol: a report of the American College of Cardiology/American Heart Association Task Force on Clinical Practice Guidelines. Circulation. 2019;139:e1082–143.
Diabetes Genetics Initiative of Broad Institute of Harvard and MIT, Lund University, and Novartis Institutes of BioMedical Research, Saxena R, Voight BF, Lyssenko V, Burtt NP, de Bakker PI, Chen H, Roix JJ, et al. Genome-wide association analysis identifies loci for type 2 diabetes and triglyceride levels. Science. 2007;316:1331–6.
Kathiresan S, Manning AK, Demissie S, D’Agostino RB, Surti A, Guiducci C, Gianniny L, Burtt NP, Melander O, Orho-Melander M, et al. A genome-wide association study for blood lipid phenotypes in the Framingham Heart Study. BMC Med Genet. 2007;8(Suppl 1):S17.
Kathiresan S, Melander O, Anevski D, Guiducci C, Burtt NP, Roos C, Hirschhorn JN, Berglund G, Hedblad B, Groop L, et al. Polymorphisms associated with cholesterol and risk of cardiovascular events. N Engl J Med. 2008;358:1240–9.
Teslovich TM, Musunuru K, Smith AV, Edmondson AC, Stylianou IM, Koseki M, Pirruccello JP, Ripatti S, Chasman DI, Willer CJ, et al. Biological, clinical and population relevance of 95 loci for blood lipids. Nature. 2010;466:707–13.
Asselbergs FW, Guo Y, van Iperen EP, Sivapalaratnam S, Tragante V, Lanktree MB, Lange LA, Almoguera B, Appelman YE, Barnard J, et al. Large-scale gene-centric meta-analysis across 32 studies identifies multiple lipid loci. Am J Hum Genet. 2012;91:823–38.
Albrechtsen A, Grarup N, Li Y, Sparso T, Tian G, Cao H, Jiang T, Kim SY, Korneliussen T, Li Q, et al. Exome sequencing-driven discovery of coding polymorphisms associated with common metabolic phenotypes. Diabetologia. 2013;56:298–310.
Tachmazidou I, Dedoussis G, Southam L, Farmaki AE, Ritchie GR, Xifara DK, Matchan A, Hatzikotoulas K, Rayner NW, Chen Y, et al. A rare functional cardioprotective APOC3 variant has risen in frequency in distinct population isolates. Nat Commun. 2013;4:2872.
Willer CJ, Schmidt EM, Sengupta S, Peloso GM, Gustafsson S, Kanoni S, Ganna A, Chen J, Buchkovich ML, Mora S, et al. Discovery and refinement of loci associated with lipid levels. Nat Genet. 2013;45:1274–83.
Holmen OL, Zhang H, Fan Y, Hovelson DH, Schmidt EM, Zhou W, Guo Y, Zhang J, Langhammer A, Lochen ML, et al. Systematic evaluation of coding variation identifies a candidate causal variant in TM6SF2 influencing total cholesterol and myocardial infarction risk. Nat Genet. 2014;46:345–51.
Peloso GM, Auer PL, Bis JC, Voorman A, Morrison AC, Stitziel NO, Brody JA, Khetarpal SA, Crosby JR, Fornage M, et al. Association of low-frequency and rare coding-sequence variants with blood lipids and coronary heart disease in 56,000 whites and blacks. Am J Hum Genet. 2014;94:223–32.
Surakka I, Horikoshi M, Magi R, Sarin AP, Mahajan A, Lagou V, Marullo L, Ferreira T, Miraglio B, Timonen S, et al. The impact of low-frequency and rare variants on lipid levels. Nat Genet. 2015;47:589–97.
Tang CS, Zhang H, Cheung CY, Xu M, Ho JC, Zhou W, Cherny SS, Zhang Y, Holmen O, Au KW, et al. Exome-wide association analysis reveals novel coding sequence variants associated with lipid traits in Chinese. Nat Commun. 2015;6:10206.
van Leeuwen EM, Karssen LC, Deelen J, Isaacs A, Medina-Gomez C, Mbarek H, Kanterakis A, Trompet S, Postmus I, Verweij N, et al. Genome of The Netherlands population-specific imputations identify an ABCA6 variant associated with cholesterol levels. Nat Commun. 2015;6:6065.
Iotchkova V, Huang J, Morris JA, Jain D, Barbieri C, Walter K, Min JL, Chen L, Astle W, Cocca M, et al. Discovery and refinement of genetic loci associated with cardiometabolic risk using dense imputation maps. Nat Genet. 2016;48:1303–12.
Liu DJ, Peloso GM, Yu H, Butterworth AS, Wang X, Mahajan A, Saleheen D, Emdin C, Alam D, Alves AC, et al. Exome-wide association study of plasma lipids in >300,000 individuals. Nat Genet. 2017;49:1758–66.
Lu X, Peloso GM, Liu DJ, Wu Y, Zhang H, Zhou W, Li J, Tang CS, Dorajoo R, Li H, et al. Exome chip meta-analysis identifies novel loci and East Asian-specific coding variants that contribute to lipid levels and coronary artery disease. Nat Genet. 2017;49:1722–30.
Hoffmann TJ, Theusch E, Haldar T, Ranatunga DK, Jorgenson E, Medina MW, Kvale MN, Kwok PY, Schaefer C, Krauss RM, et al. A large electronic-health-record-based genome-wide study of serum lipids. Nat Genet. 2018;50:401–13.
Kanai M, Akiyama M, Takahashi A, Matoba N, Momozawa Y, Ikeda M, Iwata N, Ikegawa S, Hirata M, Matsuda K, et al. Genetic analysis of quantitative traits in the Japanese population links cell types to complex human diseases. Nat Genet. 2018;50:390–400.
Klarin D, Damrauer SM, Cho K, Sun YV, Teslovich TM, Honerlaw J, Gagnon DR, DuVall SL, Li J, Peloso GM, et al. Genetics of blood lipids among ~300,000 multi-ethnic participants of the Million Veteran Program. Nat Genet. 2018;50:1514–23.