We collected the summary statistics from the largest published GWAS to maximize our power for novel discovery. PAD summary statistics were the European ancestry subjects from the recent VA Million Veteran Program analysis which consisted of 24,009 PAD cases and 150,983 PAD controls (Klarin et al., 2019). These data can be accessed from dbGAP (phs001672). CAD data were taken from CARDIoGRAMplusC4D combined with the United Kingdom BioBank (UKBB) (Van Der Harst and Verweij, 2018) and consisted of 122,733 CAD cases and 424,528 CAD controls. Data for body mass index (BMI) (meta-analysis of GIANT and UKBB; 806,834 individuals; (Yengo et al., 2018)), type 2 diabetes (T2D) (meta-analysis of consortia; 228,499 cases and 1,178,783 controls; (Vujkovic et al., 2020)), smoking initiation (smoking) (UKBB; 462,690 individuals; (Wootton et al., 2020)), and 4 lipid traits (meta-analysis of MVP and GLGC data; 723,000 participants; (Klarin et al., 2018)). Access urls for all data obtained from the public domain are provided in Supplementary Table S1.

Using the summary statistics from publicly available single-trait GWAS (Supplementary Table S1), we performed 15 N-GWAMA (Baselmans et al., 2019) multi-trait GWAS centered around PAD, CAD, and the following atherosclerotic risk factor traits: BMI, smoking, T2D, LDL-C, HDL-C, TC, and TG. Full details are provided in Supplementary Methods. Briefly, we first performed a bivariate GWAS for PAD and CAD followed by a series of trivariate GWAS combining PAD, CAD, and one of seven correlated traits that represented traditional atherosclerotic risk factors. We also performed a series of bivariate GWAS between PAD and these seven traits individually, given that a series of bivariate GWAS between CAD and most of these seven traits has already been performed (Siewert and Voight, 2018).

Each N-GWAMA multi-trait GWAS resulted in a set of independent loci represented by a sentinel SNP. We defined an independent locus as the genomic region that includes all variants within 1 megabase (Mb) of the sentinel SNP and any other variants that were in linkage diseqilibrium (LD) of r ² > 0.2 with the sentinel SNP using the 1,000 Genomes European ancestry cohort (1 kG EUR) (Purcell et al., 2007). We then applied a series of filters to remove loci that were not plausibly pleiotropic, or did not represent novel associations. To ensure there was evidence a locus was pleiotropic and that a single trait was not driving the association, we required that the sentinel SNP was at least nominally associated (p < 5 × 10⁻³) with all the individual traits involved in the multi-trait analysis. We also required that none of the variants at an independent locus were previously associated with any of the traits used in the N-GWAMA multi-trait GWAS by applying two filters. First, we required each sentinel SNP was not genome-wide significant for any of the individual traits (p > 5 × 10⁻⁸). It was also necessary that none of the SNPs at the independent locus were previously reported to be genome-wide significant for any of the individual traits involved in the multi-trait analysis in the GWAS Catalog (Buniello et al., 2019) (Figure 1). Finally, we excluded loci in the HLA region from these experiments due to the difficulty of interpreting the independent signals of these loci. Code for the pipeline is available at: https://github.com/Bellomot/Athero_NGWAMA_Multitrait_GWAS.

Given that we performed 15 multi-trait GWAS for combinations of related traits, we next implemented a multiple testing correction procedure to assess significance. Due to the high correlation between each of the multi-trait GWAS we performed, a Bonferroni correction (p < 3.3 × 10⁻⁹) for each trait combination test is conservative. Thus, we constructed a null distribution Z-score sampling strategy to estimate an α = 0.05 p-value threshold given the set of N-GWAMA multi-trait GWAS that we performed. Under the assumption that the correlation of the Z-scores across the entire genome that resulted from the N-GWAMA multi-trait GWAS are a reasonable estimate of the correlation in the multivariate null distribution of Z-scores, we can use the correlation matrix of the Z-scores from the 15 multi-trait GWAS to model the 15-dimentional multivariate standard normal distribution that is the theoretical null distribution of these results.

To get this estimate of the appropriate α = 0.05 p-value threshold, we first drew 10,000 sets of 1 million samples from a 15-dimentional multivariate normal distribution centered at the origin, and used the correlation between the Z-scores of all the SNPs that were tested across all 15 N-GWAMA multi-trait GWAS as the correlation matrix (Supplementary Table S2). We kept the most extreme Z-score from each of the 10,000 sets and then identified the 95th percentile of the most extreme Z-scores as our α = 0.05 threshold. We defined experiment-wide significance as the 95th percentile Z-score of 5.87, which corresponds to a p-value of 4.3 × 10⁻⁹.

For each multi-trait GWAS, we assessed the evidence of a shared causal variant at each significant locus by performing colocalization analysis between the trait signals using COLOC for bivariate GWAS and MOLOC for trivariate GWAS (Giambartolomei et al., 2014; Giambartolomei et al., 2018) (Supplementary Table S3). For this analysis, we applied a 500 kilobases (Kb) window (+/− 250 Kb) around the sentinel SNP. A conditional probability of colocalization is defined as the posterior probability of colocalization conditioned on the presence of a signal for each trait. For example, when using COLOC, this would be represented as the posterior probability of hypothesis four (PP4) over the sum of the posterior probabilities of hypothesis three (PP3) and PP4 (PP4/(PP3 + PP4)) (Giambartolomei et al., 2014). A probability of ≥0.8 was considered of interest. Loci that had a conditional probability of colocalization >0.5 and <0.8 were visually inspected using LocusZoom plots (Supplementary Table S4). If the LD structure suggested additional associations unlinked to the sentinel SNP, approximate conditional analysis was performed (see details on Approximate Conditional Analysis below).

We performed single-tissue colocalization analysis to prioritize candidate causal genes implicated in each N-GWAMA multi-trait GWAS using RNA-seq data obtained from the Genotype-Tissue Expression (GTEx) project (Giambartolomei et al., 2014). We identified the list of genes and tissues for which each N-GWAMA sentinel SNP was a significant single-tissue eQTL in any GTEx v8 tissue as reported in the “.signif_variant_gene_pairs.txt.gz” files available from the GTEx Portal (Supplementary Table S5). We then performed colocalization between either CAD or PAD, as determined by which trait had the most significant sentinel SNP at each locus, and each single-tissue eQTL signal (Giambartolomei et al., 2014). We selected the window of colocalization to be 500 Kb spanning the sentinel SNP. Similar to trait-trait colocalization analysis, our threshold to classify the traits as colocalized was a conditional probability of colocalization (PP4/(PP3 + PP4)) ≥ 0.8. We visually inspected LocusZoom plots for loci where colocalization analysis resulted in a conditional probability of colocalization <0.8 but >0.5 and performed approximate conditional analysis when the LD structure suggested possible allelic heterogeneity (see details on Approximate Conditional Analysis below).

For each locus that showed evidence of multiple independent signals, we performed approximate conditional analysis on variants that appeared to be associated with the trait of interest independently of the sentinel SNP (Supplementary Table S6). This analysis was necessary given that the presence of multiple associated variants in a region violates the assumptions of COLOC (Giambartolomei et al., 2014). Potential nearby association signals were identified using LocusZoom plots and the LDlink LDassoc tool between sentinel and putative distinctive signal variant quantified using LDlink (Pruim et al., 2010; Machiela and Chanock, 2018). We performed approximate conditional analysis using GCTA-COJO with 1,000 Genome Project data (European samples, n = 503) as a reference panel (Yang et al., 2011; 1000 Genomes Project Consortium et al., 2012). We conditioned the sentinel SNP on the most associated variant for each potential confounding signal identified at the locus. We then repeated the colocalization experiment on the locus using the conditional variant p-values. A full list of traits, the sentinel SNPs, and the conditioned variants for each conditional analysis are provided in the supplement (Supplementary Table S6).

We performed a colocalization analysis between the PAD signal at the PCSK6 locus and the GTEx v8 liver tissue splicing quantitative trait locus (sQTL) signal with the intron ID: 101365044:101366196:clu_14775. We also identified this intron signal in the Phenotyping Lipid traits in iPS derived hepatocytes Study (PhLiPS) hepatocyte-like cell (HLC) sQTL data by lifting over the start and stop of this intron to hg19 (101905249:101906401) (Gawronski et al., 2019). We then performed colocalization analysis between the HLC sQTL signal and the PAD signal as well as the HLC sQTL signal and the GTEx v8 liver tissue sQTL signal (Figure 2).

We first calculated the genetic correlation between PAD and CAD and the seven atherosclerosis risk factors (BMI, smoking, T2D, LDL-C, HDL-C, TC, and TG) using the summary statistics files for these GWAS using LD score regression (Supplementary Figure S1; Supplementary Table S7). We then performed 15 N-GWAMA scans centered around PAD and CAD to detect novel loci not previously reported as genome-wide significant for any of their respective traits (Supplementary Table S8). A total of 150 sentinel SNPs were multivariate genome-wide significant with all single trait p-values between <5 × 10⁻³ and >5 × 10⁻⁸ (Supplementary Table S3). Of these sentinel SNPs, 31 were nominal genome-wide significant (multi-trait p < 5 × 10⁻⁸) and met our trait-to-trait colocalization criteria, and thirteen were experiment-wide significant (multi-trait p < 4.3 × 10⁻⁹) and met our trait-to-trait colocalization criteria. If we had used the Bonferroni correction threshold, only one locus, NUP85 (p = 3.55 × 10⁻⁹), would change from experiment-wide significant to genome-wide significant. The nominal genome-wide significant sentinel SNPs represent 25 independent loci and 11 experiment-wide significant independent loci (>1 Mb from any of our other reported loci, Table 1, Supplementary Figures S2–S34). Fourteen of the nominal genome-wide significant sentinel SNPs and eight experiment-wide significant sentinel SNPs colocalized with one or more single-tissue eQTLs (Supplementary Table S5). Finally, five loci colocalized with eQTLs for genes that have been implicated in atherosclerosis by previous studies.

We noted that two of our signals that exceeded multi-trait experiment-wide significance mapped to a previously established locus for several cardiometabolic traits. That signal was tagged by rs2076211, and associated with PAD, CAD, and BMI (bivariate p = 7.6 × 10⁻¹⁰) or PAD, CAD, and LDL (bivariate p = 4 × 10⁻⁹). This variant mapped to the nearby gene PNPLA3, a well-established locus associated for non-alcoholic fatty liver disease (Speliotes et al., 2011), multiple liver enzymes measures (Yuan et al., 2008), hemotological traits (Kamatani et al., 2010), sex-hormone binding globulin levels (Ruth et al., 2020), and T2D (Mahajan et al., 2018). Our sentinel SNP that tagged both multi-trait signals was in strong LD with rs738409 (r ² = 0.73 1 kG EUR), the previously reported sentinel variant associated with these additional traits. A previous multi-trait scan for T2D and LDL reported strong association for both traits at this locus (Klimentidis et al., 2020). However, association with this locus and atherosclerotic disease (PAD, CAD) or to BMI to our knowledge has not been previously reported, but are compelling given the extensive pleiotropy for atherosclerotic causal risk factors here.

We detected a signal that exceeded multi-trait experiment-wide significance with PAD and LDL-C (bivariate p = 3.2 × 10⁻¹⁰) at the PCSK6 locus. A rare coding variant in this region has been reported to associate with LDL-C (Klimentidis et al., 2020; Sinnott-Armstrong et al., 2021), however, the coding variant (NP_002561.1:p.Thr964Met, rs34631529) and our sentinel SNP (rs1531817) are not in linkage disequilibrium (r ² = 0.0086 1 kG EUR) based on data from the 1,000 Genomes Project (Purcell et al., 2007), indicating that we detected a different signal at this locus. To further differentiate whether our signal was novel, we performed an additional conditional analysis on the coding variant rs34631529 in PAD data without any notable changes in the PCSK6 locus signal (Supplementary Figure S36). We further note that previous GWAS have found that variants at this locus are associated with inflammatory markers (Iyengar et al., 2015; Hackinger et al., 2018; Nath et al., 2019; Folkersen et al., 2020; Richardson et al., 2020; Ruotsalainen et al., 2021).

To better understand how genetic variation at the PCSK6 locus influences circulating lipid levels, we investigated the association of the bivariate sentinel SNP at this locus in the publicly available GWAS of NMR lipid subfractions: extra-small subfrations (XS), extra-large subfractions (XL), HDL, intermediate density lipoprotein (IDL), LDL, and very-low density lipoprotein (VLDL) (Kettunen et al., 2016). We found our sentinel SNP (rs1531817) had a nominal association with medium VLDL particles (β = 0.03, SE = 0.01, p = 9 × 10⁻³), total lipids in medium VLDL (β = 0.03, SE = 0.01, p = 0.02), TG in medium VLDL (β = 0.02, SE = 0.01, p = 0.03), serum TG (β = 0.02, SE = 0.01, p = 0.03), and TG in large VLDL (β = 0.02, SE = 0.01, p = 0.03).

Our sentinel SNP was a sQTL for PCSK6 in GTEx v8 liver tissue (Figure 2). To identify a potential experimental model of this splicing change, we searched for this sQTL in PhLiPS HLC summary data (PhLiPS Study) (Pashos et al., 2017). The signal at PCSK6 colocalized with an sQTL in these data as well (Figure 2), which suggests that derived hepatocyte-like cells could be a good model for further studied of the effect of this locus on LDL-C and PAD risk.

We detected a nominal genome-wide significant signal with PAD and TC (bivariate p = 2.4 × 10⁻⁸) at the SORCS3 locus rs11599236 (Table 1). This locus was previously observed to be genome-wide significant in GWAS studies of mood disorders (Howard et al., 2018; Ward et al., 2020) (Supplementary Table S4). This signal colocalized with SORCS3 mRNA expression in pituitary tissue (Supplementary Figure S34). The opposite direction of effect was noted for both traits and the gene-tissue pair: decreased SORCS3 associated with increased PAD and TC (β = 0.31, SE = 0.05, p = 5.5 × 10⁻⁸).

We detected a trivariate GWAS signal with PAD, CAD, and smoking (trivariate p = 1.9 × 10⁻⁹) at the HMBS locus rs1006195 (Table 1). This variant was genome-wide significant in previous GWAS studies for several cardiometabolic traits including Apolipoprotein A1 levels, waist-hip ratio, BMI, fat mass percentage, HDL-C, and T2D (Turcot et al., 2018; Van Der Harst and Verweij, 2018; Pulit et al., 2019; Richardson et al., 2020; Zhu et al., 2020). This pleiotropic signal colocalized with HMBS and VPS11 mRNA expression in several tissues (Supplementary Figure S14; Supplementary Table S5). HMBS (β = 0.28, SE = 0.03, p = 4.1 × 10⁻²³) demonstrated the same direction of effect with PAD, CAD, and smoking, indicating that increased HMBS expression is associated with increased PAD, CAD, and smoking risk. VPS11 (β = −0.25, SE = 0.03, p = 6.0 × 10⁻¹⁵) demonstrated the opposite direction of effect in all tissue except skeletal muscle and the left ventricle of the heart, meaning increased VPS11 expression in skeletal muscle and left ventricle is associated with increased PAD, CAD, and smoking risk.

We also detected a nominal genome-wide significant signal with PAD and CAD (bivariate p = 7.8 × 10⁻⁹) at the VDAC2 locus rs7088974 (Table 1). Variants at this locus have been found to be associated with BMI in previous GWAS(Pulit et al., 2019) (Supplementary Table S4). Although variants near this locus have also been associated with smoking behavior, our results suggest that the locus we detected is independent of smoking behavior (Liu et al., 2019) (Supplementary Figure S35). This signal colocalized with VDAC2 mRNA expression in multiple vascular tissues relevant to atherosclerosis, including aorta and tibial artery (Supplementary Table S5). The direction of effect in all tissue was opposite to the direction of effect for PAD and CAD: the allele associated with increased VDAC2 expression (EA = C, EAF = 0.57) is associated with decreased PAD and CAD (β = −0.14, SE = 0.02, p = 7.7 × 10⁻⁹).

Finally, we detected a trivariate signal with PAD, CAD, and HDL-C (trivariate p = 3.4 × 10⁻¹¹) at the LRCH1 locus rs9526214 (Table 1). This locus had evidence of allelic heterogeneity when we reviewed the regional association plots, which led us to perform approximate conditional analyses on the pleiotropic signal sentinel SNP rs9316223 and the resulting conditional probability of colocalization met our criteria (Supplementary Table S5). This locus has been found to be genome-wide significant in previous GWAS studies for platelets, systolic blood pressure, and stroke (Evangelou et al., 2018; Malik et al., 2018) (Supplementary Table S4). This signal colocalized with LRCH1 mRNA expression in tibial artery, whole blood, and other tissues (Supplementary Figure S12; Supplementary Table S5). The opposite direction of effect was noted for all three traits and the gene tissue pair: the allele associated with increased LRCH1 was also associated with decreased PAD, CAD, and HDL-C (β = −0.09, SE = 0.02, p = 1.3 × 10⁻⁸).

We identified a signal at the PCSK6 locus that has a bivariate association with PAD and LDL-C and provide strong evidence in support of PCSK6 as the causal gene at the PCSK6 locus. PCSK6 is a calcium-dependent serine endoprotease that cleaves proteins to active and inactive forms depending on the target protein (Kiefer et al., 1991).

There is convincing experimental evidence to suggest that PCSK6 directly influences plaque development and plaque stability. Smooth muscle cell migration in injured arteries is facilitated by cytokine induced PCSK6 expression that activates matrix metalloproteinases (MMP14/MMP2) (Perisic et al., 2013). This smooth muscle cell mechanism may explain the association of PCSK6 with carotid intima-media thickness in a candidate gene study (Rykaczewska et al., 2020). The sentinel SNP associated with maximum progression of carotid intima-media thickness was the same as the sentinel SNP identified in our bivariate scan between LDL-C and PAD. There is also evidence that PCSK6 activates MMP9, which enhances degredation of the extra cellular matrix and thus promotes plaque instability (Li et al., 2020; Testa et al., 2021).

Our data suggest PCSK6 also influences lipid metabolism, a known upstream cause of atherosclerotic progression. This is in agreement with the known role of PCSK6 in lipid metabolism, where it cleaves and inactivates endothelial lipase (EL) and lipoprotein lipase (LPL) (Jin et al., 2005), which can lead to hyperlipidemia (Choi and Korstanje, 2013). This finding is further supported by the recently reported rare coding variant in PCSK6 associated with decreased LDL-C (Klimentidis et al., 2020; Sinnott-Armstrong et al., 2021).

It remains to be determined how changes in PCSK6 activity result in altered lipid metabolism and risk of PAD. It is possible that changes in PCSK6 activity alter lipid metabolism and therefore have downstream effects on PAD or altered PCSK6 activity may effect LDL-C and PAD separately through independent mechanisms. In all likelihood, it is a combination of LDL-C dependent and independent mechanisms that link PCSK6 to PAD.

From a drug targeting standpoint, we would expect that inhibiting PCSK6 reduces LDL-C levels and PAD risk. There are several non-FDA approved, non-specific PCSK6 inhibitors that influence lipid metabolism: alpha1-antitrypsin Portland (alpha1-PDX) (Gauster et al., 2005), profurin (Jin et al., 2007), Pf-pep (Byun et al., 2010), and dicoumarols (DC), specifically DC2 (Komiyama et al., 2009).

We identified a novel pleiotropic association at the SORCS3 locus between PAD and TC. Colocalization experiments support SORCS3 as the causal gene at this locus. SORCS3 is a type I transmembrane protein that is a member of the Vps10p receptor family (Hermey, 2009). The post-synaptic sorting receptor SORCS3 is highly expressed in the hippocampus and binds NGF and PDGF-BB to modulate several signal transduction pathways involved in neuronal activity (Hermey, 2009; Christiansen et al., 2017). A recent knockout mouse study has shown that even loss of one Sorcs3 allele decreases lipid metabolism as a source of energy and increases adiposity (Subkhangulova et al., 2018). The proposed mechanism for this occurrence is that loss of SORCS3 expressed in agouti-related peptide (AGRP) neurons leads to unchecked production of circulating AGRP (Henry et al., 2015). ARGP is well described to cause enhanced food intake, decreased locomotor activity, decreased use of lipids as metabolic fuel, and overall increased adiposity (Maier et al., 2018; Beutler et al., 2020). One limitation of GTEx data is a very small sample size of hypothalamus data, our data show that reduced expression of SORCS3 in the pituitary gland is related to high levels of TC and an increased overall risk for PAD, potentially through increased circulating AGRP diverting energy metabolism away from lipid sources.

We acknowledge there are several limitations to this study. First, there is sample overlap between several of our single trait summary statistics files. The N-GWAMA method attempts to account for this, but if the correction for the overlap was insufficient this could inflate our false discovery rate. Second, in using the null Z-score resampling strategy to estimate the appropriate multiple testing correction, we assume that the correlation of the Z-scores across the entire genome that resulted from the N-GWAMA multi-trait GWAS are a reasonable approximation of the correlation between the multivariate null distribution of Z-scores. This assumption could be broken due to highly associated SNPs leading to an altered correlation matrix compared to the true null distribution and the Z-scores having an inflated variance from confounders that may not have been corrected for by LDSC. Finally, we selected nine atherosclerotic and cardiometabolic traits based on conventional relationships with atherosclerosis; however, there are likely multifactorial and multidirectional relationships within this group of traits. It is possible that some of the novel loci represent the interaction between traits instead of the intended representation of atherosclerosis as we have interpreted it. It is also possible that including other sets of cardiometabolic risk factors may identify additional novel loci.