Genetically Predicted Glucose-Dependent Insulinotropic Polypeptide (GIP) Levels and Cardiovascular Disease Risk Are Driven by Distinct Causal Variants in the GIPR Region.

There is considerable interest in GIPR agonism to enhance the insulinotropic and extrapancreatic effects of GIP, thereby improving glycemic and weight control in type 2 diabetes (T2D) and obesity. Recent genetic epidemiological evidence has implicated higher GIPR-mediated GIP levels in raising coronary artery disease (CAD) risk, a potential safety concern for GIPR agonism. We therefore aimed to quantitatively assess whether the association between higher GIPR-mediated fasting GIP levels and CAD risk is mediated via GIPR or is instead the result of linkage disequilibrium (LD) confounding between variants at the GIPR locus. Using Bayesian multitrait colocalization, we identified a GIPR missense variant, rs1800437 (G allele; E354), as the putatively causal variant shared among fasting GIP levels, glycemic traits, and adiposity-related traits (posterior probability for colocalization [PPcoloc] > 0.97; PP explained by the candidate variant [PPexplained] = 1) that was independent from a cluster of CAD and lipid traits driven by a known missense variant in APOE (rs7412; distance to E354 ∼770 Kb; R 2 with E354 = 0.004; PPcoloc > 0.99; PPexplained = 1). Further, conditioning the association between E354 and CAD on the residual LD with rs7412, we observed slight attenuation in association, but it remained significant (odds ratio [OR] per copy of E354 after adjustment 1.03; 95% CI 1.02, 1.04; P = 0.003). Instead, E354's association with CAD was completely attenuated when conditioning on an additional established CAD signal, rs1964272 (R 2 with E354 = 0.27), an intronic variant in SNRPD2 (OR for E354 after adjustment for rs1964272: 1.01; 95% CI 0.99, 1.03; P = 0.06). We demonstrate that associations with GIP and anthropometric and glycemic traits are driven by genetic signals distinct from those driving CAD and lipid traits in the GIPR region and that higher E354-mediated fasting GIP levels are not associated with CAD risk. These findings provide evidence that the inclusion of GIPR agonism in dual GIPR/GLP1R agonists could potentiate the protective effect of GLP-1 agonists on diabetes without undue CAD risk, an aspect that has yet to be assessed in clinical trials.

Introduction

The incretin hormones glucose-dependent insulinotropic polypeptide (GIP) and glucagon-like peptide 1 (GLP-1) are well-known for their insulinotropic activity (1,2), which is diminished in type 2 diabetes (T2D) (3–6). This has prompted significant therapeutic interest in the agonism of their respective receptors, GIPR and GLP1R, to enhance their insulinotropic and extrapancreatic effects (7,8). Moreover, preclinical and clinical data demonstrate that dual agonism of GIPR and GLP1R delivers superior glycemic and weight control efficacy compared with selective GLP1R agonism (9–12). Clinical proof for the superiority of tirzepatide, a dual GIPR/GLP1R agonist, versus GLP1R agonism was established in a 6-month dose-range-finding phase 2b trial in subjects with T2D (11). Post hoc analysis reported a beneficial effect on cardiovascular risk biomarkers compared with the blinded GLP1R agonist included in the trial (13,14).

Little direct preclinical experimental evidence exists for GIPR agonism contributing to cardiovascular disease (CVD) risk (15,16). GIP exhibits antiatherogenic effects on vascular endothelial cells (17–20) with the exception that it has been reported to stimulate expression of osteopontin in the vasculature in an endothelin-1–dependent manner (21). Additionally, GIP exerts anti-inflammatory effects on monocytes/macrophages (17,22). These in vitro findings are reflected by cardioprotective GIP pharmacology in mouse models of atherosclerosis irrespective of their diabetes condition (17,22,23). Further, GIP infusion or overexpression is protective in mouse models of restenosis and cardiac remodeling (17,24). While germline or cardiomyocyte-selective knockout of GIPR protected against ischemic injury, GIP itself was not deleterious (25). Further, cardiac selective knockout of the GIPR was not protective in experimental models of heart failure (25). In contrast with these preclinical experimental findings, recent evidence suggests that fasting GIP levels are associated with increased carotid intimal thickening (26). In addition, evidence from a recent meta-analysis (27) of two large population-based cohort studies suggests that higher fasting but not postchallenge GIP levels were associated with increased risk of CVD mortality (hazard ratio 1.30; 95% CI 1.11, 1.52; P = 0.001). GLP-1 was not associated with CVD mortality, consistent with clinical trial data (28–31) and genetic evidence (32) highlighting the beneficial effects of GLP1R agonism.

Genetic evidence from two-sample Mendelian randomization (2SMR) has reinforced suggestions that higher GIP levels raise CVD risk (27). A missense variant in GIPR, rs1800437 (E354Q), encoding a substitution of glutamic acid for glutamine at position 354 of the GIPR protein, was used as an instrumental variable for fasting GIP levels (27). The 354Q allele has been reported to reduce GIPR signaling by increasing receptor desensitization and downregulation (33). This variant has previously been associated with higher 2-h glucose (34), BMI (35), and fasting and 2-h GIP levels (36). In line with a predicted causal direction from fasting GIP levels to coronary artery disease (CAD) risk, estimates in the reverse direction showed no significant effect of CAD on fasting GIP levels (27). These estimates should be interpreted with caution, however, as 1) they represent the association of a single variant with CAD risk and do not model the effects of other variants in the region, which may dampen or modulate this effect, and 2) they do not take into account that the association between E354 and CAD may be entirely synthetic due to linkage disequilibrium (LD) between this variant and the true CAD causal variant.

Considering the pharmacological interest in modulating this pathway as a potential T2D therapeutic, increases in CVD risk would represent a major concern regarding the safety and continued development of these therapies. We aimed to quantitatively assess whether the association between higher GIPR-mediated fasting GIP levels and CAD risk is mediated via GIPR or the result of LD between variants in GIPR and other variants in the region. Using 2SMR, we aimed to quantify the association of higher fasting GIP levels with CAD and other metabolically relevant traits, including ∼6,000 omics biomarkers, using E354 as an instrumental variable. Next, using Bayesian colocalization, we aimed to partition the traits associated with E354 into distinct clusters driven by shared independent variants. Finally, using conditional analysis we aimed to assess whether any of these associations are confounded by LD between E354 and other variants in the GIPR region.

Research Design and Methods

Study Design

Three sets of genetic analyses were used to investigate the relationship between higher GIPR-mediated fasting GIP levels and CVD risk. Firstly, using univariate 2SMR, we explored the association of higher fasting GIP levels with CAD and 23 different cardiometabolic diseases, along with anthropometric, glycemic, and lipid traits and ∼6,000 omics biomarkers from both in-house and publicly available data, with E354 as a proxy (Supplementary Table 1). Next, Bayesian multitrait colocalization was used to partition the traits associated with E354 into distinct clusters driven by shared causal variants. Finally, conditional analyses were used to assess whether any of the associations with E354 are confounded by LD between E354 and other variants in the GIPR region, implying that their associations are mediated via not GIPR but, rather, other genes in the region.

Study Participants

European Prospective Investigation into Cancer and Nutrition (EPIC)-Norfolk (37) (Supplementary Table 2) is a population-based prospective cohort of individuals aged between 40 and 79 years living in Norfolk (a county of the U.K.) at the time of recruitment from primary care outpatient clinics in the city of Norwich and surrounding areas. EPIC-Norfolk (37) consists of two subcohorts, a T2D case-cohort and a quasi-random selection of participants from the larger EPIC (38,39) study. The study was approved by the Norfolk Research Ethics Committee (reference no. 05/Q0101/191), and all participants gave written consent before entering the study.

Fenland (40) (Supplementary Table 2) is a population-based cohort study of individuals without diabetes who were born between the years of 1950 and 1975 and recruited through population-based general practice registers in Cambridge, Ely, and Wisbech (Cambridgeshire County, U.K.). Ethics approval for the study was given by the Cambridge local ethics committee (reference no. 04/Q0108/19), and all participants gave written consent prior to entering the study.

UK Biobank (41) (Supplementary Table 2) is a population-based cohort study of individuals recruited from 22 rural and urban recruitment centers in the U.K. European ancestry participants with available genome-wide genotyping and phenotypic data were included in this study. Ethics approval for the UK Biobank study was given by the North West - Haydock Research Ethics Committee (16/NW/0274). This research was conducted using application 44448. Participants gave electronic consent for use of their anonymized data and samples for health-related research, to be recontacted for further substudies, and for access to their health-related records.

Genotyping and Imputation

Genome-wide genotyping in the Fenland cohort was performed in three subcohorts with use of the Affymetrix Genome-Wide Human SNP Array 5.0, the Affymetrix UK Biobank Axiom Array, or the Illumina CoreExome-24 v1 BeadChip, with imputation to the Haplotype Reference Consortium v1.1 (42), the 1000 Genomes Project (43), and the UK10K (44) reference panels. Samples from EPIC-Norfolk and UK Biobank were genotyped with the Affymetrix UK Biobank Axiom Array and imputed to the same reference panels.

Profiling of the Plasma Proteome

Fasting EDTA plasma samples from 12,084 participants from the Fenland (40) study were subjected to proteomic profiling by SomaLogic (Boulder, CO) using an aptamer-based technology (somascan v4). The relative abundances of 4,775 human proteins were measured using 4,979 SOMAmers (45). For accounting for within-run hybridization variability, control probes were used to generate a scaling factor for each sample. Differences in total signal between samples as a result of variation in overall protein concentration or technical variability such as reagent concentration, pipetting, or assay timing were accounted for using the ratio between each SOMAmer measured value and a reference value. The median of these ratios was computed for each dilution set (40%, 1%, and 0.005%) and applied to each dilution set. Samples were removed if they failed SomaLogic quality control measures or did not meet the acceptance criteria of between 0.25 and 4.00 for all scaling factors. A total of 10,078 samples had available genotype data and were used in this study. Aptamer target annotations and mapping to UniProt accession numbers as well as gene identifiers were provided by SomaLogic.

Plasma Metabolomic Profiling

Within EPIC-Norfolk (described previously) (37), the levels of up to 1,504 metabolites were measured in three batches using the Metabolon DiscoveryHD4 platform (46) (Metabolon, Durham, NC), in citrate plasma samples collected at baseline. Measurements were made in ∼12,000 samples, in two sets of ∼6,000 quasi–randomly selected samples, which were preceded by measurements in an incident T2D case-cohort (N = 1,503; 857 in the subcohort).

Briefly, raw data were extracted and peaks were identified and assessed for quality by Metabolon. Metabolite identification was done by comparing measures with a curated library containing the retention time, mass-to-charge ratio, and chromatographic data of known metabolites. Each metabolite was then quantified with an area under the curve method and the data were normalized to correct for instrument tuning variations across run days. For data normalization for each run day the median value for each metabolite was set to 1, normalizing each measurement proportionately. Metabolite annotations and pathway classifications are as reported by Metabolon.

Statistical Analysis

Genome-Wide Association Study of Plasma Proteins and Pairwise Colocalization of GIP Levels With Cardiometabolic Traits

Genome-wide association study (GWAS) was performed as described in Supplementary Table 3. Two SOMAmers targeted circulating GIP, namely, 16292-288 and 5755-29. SOMAmer 16292-288 was selected against amino acids 1–93 of the precursor protein (UniProt identifier P09681), whereas 5755-29 targeted amino acids 22–153. SOMAmers are relative measures of GIP abundance; therefore, to ascertain whether the underlying genetics at GIPR were comparable with previous results (36), we performed pairwise genetic colocalization analyses between GIP measures and cardiometabolic traits.

T2D, coronary heart disease, BMI, and 2-h glucose adjusted for BMI and LDL were included as cardiometabolic traits of interest (Supplementary Table 1). Summary statistics from a GWAS of 2-h glucose adjusted for BMI in Fenland (Supplementary Table 3) were preferred over those from previous efforts (34), due to denser variant coverage. Using GWAS summary statistics for each trait, the 1-Mb regions either side of E354 (chromosome 19: 45181392–47181392) were extracted. Insertions and deletions as well as any variants with a standard error of 0 were removed. Effect estimates were aligned to the GIP-raising alleles. Pairwise colocalization was conducted using the COLOC (47) R package. Priors, p1 and p2, the prior probabilities that a variant is associated with either trait, were set to 1 × 10−4, and p12, the probability that a single variant is associated with both traits, was set to 1 × 10−5. T2D and coronary heart disease were treated as case-control traits and all other traits as quantitative. Posterior probabilities for colocalization (PPcoloc) were considered significant if they met the following criteria: (H4 + H3 ≥0.9 and H4 / H3 ≥3), where H3 is the PP for two distinct genetic signals and H4 the PP for a shared genetic signal.

GWAS of Plasma Metabolites

GWAS was performed in two sets, for all metabolites present in at least 100 individuals in both sets. The first set consisted of up to 5,841 individuals from both the subcohort of the T2D case cohort and the first batch of quasi–randomly selected samples. The second set consisted of up to 5,698 individuals from the second batch of quasi–randomly selected samples. GWAS was performed as described in Supplementary Table 3.

Association Between E354 and Cardiometabolic and Molecular Traits

This work leveraged regional GWAS summary statistics from in-house studies and data from published studies in the 1-Mb regions either side of E354. Details on all included phenotypes can be found in Supplementary Table 1. GWAS for phenotypes derived in-house were performed as described in Supplementary Table 3. Only self-reported White European participants were included for all outcomes, except for plasma metabolite measures in EPIC-Norfolk (37), where all participants were included. However, participants in EPIC-Norfolk (37) overwhelmingly self-reported as White European.

We performed univariate 2SMR using the Wald ratio method (48) to estimate the potential causal effect of fasting GIP levels on various traits (Supplementary Table 1). Genetically predicted fasting GIP levels were used as the exposure with E354 as the instrumental variable (Human Genome Organisation [HUGO] gene: GIPR; National Center for Biotechnology Information [NCBI] transcript NM_00016 4.4 c.1060G>C; protein change, E354Q; E345 variant is encoded by the G allele). All summary statistics were aligned to the fasting GIP raising allele (G) of E354. Bonferroni-corrected significance thresholds were used to ascertain statistical significance of E354 across all outcomes.

Partial Correlations Between X-12283 and Known Metabolites

To estimate the metabolite class and putative functional pathway of X-12283, we estimated partial correlations between X-12283 levels and the levels of other metabolites measured in 11,966 participants from EPIC-Norfolk.

First, missing metabolite measures were imputed within each measurement set with use of multivariate imputation by chained equations (MICE) (49) with the R package mice v3.6.0. To ensure accurate imputation, we only considered the 883 metabolites with <50% missingness within both measurement sets. Imputation was repeated a total of 20 times, generating 20 sets of fully imputed results. Following imputation, measures were standardized (mean = 0, SD = 1). For each imputation, partial correlations between metabolite pairs were calculated with the R package GeneNet v1.2.14. Partial correlation estimates were transformed with Fisher Z transformation and the R package psych v1.9.12.31, and then pooled across the 20 imputations for each measurement set, with use of Rubin’s rules (50). Estimates for the two measurement sets were then meta-analyzed, using a fixed-effects, inverse variance–weighted method in the R package meta v4.12–0, and finally back-transformed to correlation estimates. P values were calculated with the Fisher transformed partial correlations.

Partial correlation estimates with absolute values of >0.1 were then used to draw a Gaussian graphical model in Cytoscape v3.2.1. Partial correlations were considered significant at a Bonferroni significance threshold of P ≤ 1.28 × 10−7, accounting for the 389,403 metabolite pairs tested.

Multitrait Colocalization Across Cardiometabolic Traits

Multitrait colocalization (HyPrColoc) (51) was used at the GIPR locus to 1) identify cardiometabolic traits that share a common causal variant and 2) partition clusters of cardiometabolic traits driven by distinct causal variants. HyPrColoc was run using the default variant-specific prior configuration; priors 1 and 2 were set at 1 × 10−4 and 0.02, respectively; and regional and alignment thresholds of 0.5 were used (51).

Variants were extracted and excluded from GWAS summary statistics for 26 cardiometabolic traits of interest as in the pairwise colocalizations above, and all variants in perfect LD (R2 = 1) with E354 were removed. The GIP measures considered were fasting GIP as measured by SOMAmers 16292-288 and 5755-29, as well as fasting and 2-h GIP measures from the Malmö Diet and Cancer (MDC) subcohort of Almgren et al. (36) Both the MDC and Prevalence, Prediction and Prevention of diabetes (PPP)-Botnia Study cohorts were genotyped with exome-wide arrays, thereby limiting the number of variants included in the analysis in considering variants present across all traits. MDC measures were preferred to those from either the PPP-Botnia Study subcohort or the meta-analysis of the two subcohorts due to denser variant coverage, despite the PPP-Botnia Study having a larger sample size. The anthropometric traits adjusted and unadjusted for BMI (where applicable) were BMI, waist-to-hip ratio (WHR), and hip and waist circumferences. T2D and CAD were included as disease outcomes. Glycemic measures included nonfasting glucose, HbA1c, 2-h glucose adjusted for BMI, fasting glucose adjusted for BMI, and fasting insulin adjusted for BMI. GWAS summary statistics from Fenland were used for fasting and 2-h glucose as well as fasting insulin. Finally, lipid traits included LDL, HDL, total cholesterol, triglycerides, lipoprotein A, apolipoprotein (apo)A1, and apoB.

To assess sensitivity in the number and size of clusters identified, increasingly stringent prior and threshold configurations were used. Prior 2 values of 0.02, 0.01, and 0.001, and threshold values of 0.5, 0.6, 0.7, 0.8, and 0.9, were considered. T2D and CAD were considered as binary case-control traits, and all others were considered quantitative. To estimate the posterior probability (PP) that the candidate variant is the causal variant (PPcausal), we multiplied the PPcoloc by the PP explained by the candidate variant (PPexplained). Trait clusters were reported at the recommended (51) thresholds of prior 2 = 0.02, regional and alignment thresholds = 0.9.

To account for low variant coverage in the MDC cohort, we ran a secondary analysis using the same populations, configuration, and sensitivity assessments as above, while we excluded the GIP traits measured in MDC.

Finally, heat maps based on similarity matrices estimating how often trait pairs were clustered together across all algorithm parameter choices were drawn. In addition, regional association plots were drawn for each cluster with the gassocplot R package and LD data from EPIC-Norfolk. All data analysis was performed with R version 3.6.3.

Conditional Analysis at the GIPR Locus

To determine whether the association between E354 and CAD was due to LD between E354 and other CAD lead variants in the GIPR region, we performed conditional analysis using GCTA (52) v1.93.1. Using full GWAS summary statistics for CAD (53) on chromosome 19, we implemented a stepwise selection to identify independent variants associated with CAD. Selection was performed with a threshold of P < 1 × 10−5, a threshold for collinearity between variants of 0.05, and a minor allele frequency threshold of 1%. An LD reference panel from EPIC-Norfolk was used. The association between E354 and CAD was then conditioned on each independent variant to estimate whether the association was attenuated, implying that the association was due to the residual LD between E354 and an independent variant. This was repeated for all traits associated with E354. If E354 (or a proxy variant in complete LD with E354) was identified as one of the independent variants, conditional analysis was not performed. Following this, regional association plots were generated using LocusZoom v1.2. To determine whether other variants previously found to be associated with fasting GIP levels (36) were associated with CAD, we extracted their estimates from the CAD summary statistics (53).

Data and Resource Availability

The data sets analyzed during the current study are publicly available, and links are provided in Supplementary Table 1. EPIC-Norfolk and Fenland data are available upon reasonable request via the study websites (https://www.mrc-epid.cam.ac.uk/research/studies/epic-norfolk/ and https://www.mrc-epid.cam.ac.uk/research/studies/fenland/in formation-for-researchers/). GIP measures from Almgren et al. (36) are available from the relevant corresponding author upon reasonable request. All data from UK Biobank are available to approved users upon application. No applicable resources were generated or analyzed during the current study.

Results

Characterization of a Missense Variant E354 (rs1800437) in GIPR

Among the cardiometabolic disease outcomes examined, higher E354-predicted fasting GIP levels were associated with lower T2D risk (odds ratio [OR] per copy of E354, 0.97; 95% CI 0.96, 0.99; P = 3 × 10−4) (Fig. 1), an effect that strengthened following BMI adjustment (0.93; 95% CI 0.91, 0.95; P = 3 × 10−14). In line with this, lower 2-h glucose levels were observed (2-h glucose in mmol/L per copy of E354, −0.09; 95% CI −0.11, −0.07; P = 2 × 10−15) (Fig. 1). Additionally, HbA1c levels were shown to be 0.01 SD units lower per copy of E354. E354 showed a weak positive association with nonfasted glucose levels. As this phenotype captures wide-ranging physiological responses in both the fasted and postprandial state, deconvoluting this association requires further investigation. E354 was associated with higher CAD risk (OR per copy of E354, 1.03; 95% CI 1.02, 1.05; P = 2 × 10−6) (Fig. 1) and higher levels of several lipid risk factors but lower triglyceride levels (Fig. 1). E354 was not significantly associated with other CVD subtypes in UK Biobank (Supplementary Fig. 1).

Figure 1

Associations between E354 (rs1800437) and cardiometabolic disease end points, glycemic traits, cardiovascular risk factors and lipids, anthropometric traits, and biomarkers estimated with 2SMR. A: Associations with cardiometabolic disease end points are shown in blue and are represented as ORs (95% CI) for each disease per copy of rs1800437. B: Associations with glycemic traits are shown in orange and cardiovascular and lipid traits in green, and anthropometric traits and biomarkers are shown in yellow and purple, respectively. Estimates are represented as β (95% CI) for each outcome per copy of rs1800437. All traits are in SD units aside from fasting and 2-h glucose, which are in mmol/L; fasting insulin, in log (pmol/L); and HbA1c, in mmol/mol. Fold change insulin represents the fold change in insulin levels between fasting and 2-h measures. A Bonferroni significance threshold of P ≤ 0.001 was used, accounting for the number of traits tested. *HbA1c estimates are in SD units per copy of E354. The corresponding clinical units, % (mmol/mol), are as follows: −2.15% (95% CI −2.15, −2.14) and −0.07 mmol/mol (95% CI −0.07, −0.06). adj., adjustment.

Each copy of E354 was associated with 0.03 SD higher BMI (95% CI 0.03, 0.04; P = 3 × 10−59) (Fig. 1). Similar associations were observed between E354 and higher regional anthropometric measures from bio-impedance data (Supplementary Fig. 2) as well as hip and waist circumferences and waist-to-hip ratio. In addition, significant associations were found with both higher lean and fat mass from a large GWAS based on bio-impedance data (Supplementary Fig. 2).

Of the 19 biomarkers investigated, E354 was significantly associated with lower levels of only two, namely, albumin and creatinine (albumin β in SD units per copy of E354, −0.01; 95% CI −0.02, −0.01; P = 6 × 10−6; creatinine −0.02; 95% CI −0.02, −0.01; P = 1 × 10−11) (Fig. 1).

Next, we estimated the association of E354 with the fasting levels of 4,979 human proteins from the SomaScan v4 assay. Significant associations with the levels of three proteins were found (Supplementary Fig. 3), one of these being 0.08 SD higher fasting GIP levels (95% CI 0.05, 0.11; P = 4 × 10−6) as measured according to SOMAmer 16292-288. Interestingly, in our analysis we did not find a significant association between the other GIP SOMAmer, 5755-29, and E354. Lower levels of secretoglobin family 3A member 1 (SCGB3A1) and glutaminyl-peptide cyclotransferase-like protein (QPCTL) were also found to be associated with E354. In contrast with a previous report (21), no association between E354 and osteopontin was found.

Lower levels of an unidentified metabolite, X-12283 (β in SD units per copy of E354, −0.08; 95% CI −0.12, −0.05; P = 2 × 10−5) (Supplementary Fig. 4), analyzed in 8,278 participants, were found to be significantly associated with E354. A total of 11 metabolites were significantly correlated with X-12283; of these, 6 had a partial correlation estimate with X-12283 with absolute values >0.1 (Supplementary Fig. 5). In addition to significant correlations with unknown metabolites, X-12283 was most significantly correlated with indolepropionate (correlation estimate = 0.21; P = 1 × 10−45) (Supplementary Fig. 5).

Multitrait Colocalization Across Cardiometabolic Traits at GIPR

A total of 418 genetic variants were included in the main analysis, which was limited due to the inclusion of fasting and 2-h GIP measures from MDC (36), whereas 4,996 were included in the secondary analysis (Table 1). Using the recommended prior and threshold configuration, we identified five distinct trait clusters, three of which were shared by both analyses (Table 1). Cluster similarity across all prior and threshold permutations for the two analyses is summarized in heat maps (Fig. 2). Results for all permutations for both analyses can be found in Supplementary Tables 4 and 5, respectively.

Table 1

Clusters of colocalized traits identified by the main and secondary analyses at recommended settings

		Main analysis					Secondary analysis
Locus	Candidate LD (R2)§	Colocalized traits	PPcoloc*	Candidate variant	PPexplained	N variants	Colocalized traits	PPcoloc*	Candidate variant	PPexplained	N variants
GIPR	1	LDL, CAD, total cholesterol, lipoprotein A, apoB	1	rs7412	1	418	LDL, CAD, total cholesterol, lipoprotein A, apoB	1	rs7412	1	4,996
GIPR	0.69	HbA1c, apoA1, WHRadjBMI, waist circumference adjBMI, WHR	0.99	rs4420638	1	418	HbA1c, HDL, apoA1, WHRadjBMI, waist circumference adjBMI, WHR	0.97	rs429358	1	4,996
GIPR	1	GIP SOMAmer 16292-288, fasting GIP, 2-h GIP, BMI, glucose, hip circumference, waist circumference, 2-h glucose adjBMI	0.97	rs1800437	1	418	SOMAmer 16292-288, hip circumference, 2-h glucose adjBMI	0.91	rs1800437	0.68	4,996
GIPR	NA						BMI, waist circumference	1	rs1800437	1	4,996
GIPR	NA	T2D, T2D adjBMI	0.98	rs8108269	0.99	418
GIPR	NA						Triglycerides, hip circumference adjBMI	0.98	rs5117	0.93	4,996

Blank rows for either analysis indicate a cluster not identified in the respective analysis. adj, adjusted for; GIPR, glucose-dependent insulinotropic polypeptide receptor; glucose, nonfasting glucose; NA, not applicable; variants, single nucleotide polymorphisms.

Trait clusters are reported at the recommended thresholds for HyPrColoc: prior 2 = 0.02; regional and alignment thresholds = 0.9.

The LD in R2 between the candidate variants for the main and secondary analyses, respectively.

Figure 2

Similarity heat map for each cluster at the GIPR locus across prior and threshold permutations. Traits that were estimated to colocalize are clustered together. Darker colors represent traits that were estimated to colocalize more often across prior and threshold permutations (prior 2: 0.02, 0.01, and 0.001; thresholds: 0.5, 0.6, 0.7, 0.8, and 0.9). A: Main analysis. B: Secondary analysis. 2hr, 2-hour; adj, adjusted for; circumf., circumference; Glucose, nonfasting glucose; HC, hip circumference; WC, waist circumference.

Of the clusters identified, two distinct clusters were of interest. The first, driven by rs7412, a missense variant in the apoE gene (APOE), contained CAD and lipid traits—many of which are established CVD risk factors. Both PPcoloc and PPcausal were estimated to be 1 in the two analyses, demonstrating robust evidence for colocalization (Table 1 and Supplementary Fig. 6). This robustness is further emphasized as the same cluster of traits was identified in using more stringent prior configurations (Fig. 2 and Supplementary Tables 4 and 5). A second cluster of GIP, anthropometric, and glycemic traits was driven by rs1800437 (E354) (Table 1 and Supplementary Fig. 7). The PPcoloc for both analyses showed robust evidence for colocalization (main analysis, PPcoloc = 0.97, PPexplained = 1, PPcausal = 0.97; secondary analysis, PPcoloc = 0.91, PPexplained = 0.68, PPcausal = 0.62). A second cluster of BMI and waist circumference driven by E354 was observed in the secondary analysis (Table 1). Sensitivity analyses showed that this split was an artifact of the branch and bound clustering algorithm in HyPrColoc and the single causal variant assumption (Supplementary Fig. 7). Removal of the clustering algorithm showed that BMI and waist circumference were part of the larger cluster of GIP, anthropometric, and glycemic traits driven by E354 (PPcoloc = 0.95, PPexplained = 1, PPcausal = 0.95).

Critically, these results replicate our findings using pairwise-trait colocalization at this locus, showing that fasting GIP levels and CVD risk are driven by independent variants (R2 between E354 and rs7412 = 0.004) (Table 1, Supplementary Figs 6–8, and Fig. 2). Additionally, both colocalization analyses demonstrate that the underlying genetics at GIPR are comparable between GIP levels measured by SOMAmer 16292-288 and the ELISA of previous analyses (36). Together these results robustly demonstrate that the GIP-raising and CVD risk–increasing effects at this locus are distinct (Supplementary Tables 4 and 5).

The traits of a third cluster, including a mixture of glycemic and anthropometric traits and apoA1 levels, were estimated to colocalize at rs4420638, which was in LD with rs429358 (R2 = 0.69), a missense variant in APOE identified as the candidate variant in the secondary analysis (R2 with E354 = 0.001). In the secondary analysis, HDL was also included as part of the cluster. As the secondary analysis included more variants and therefore had greater genomic context, rs429358 is likely to be the candidate variant at which these traits colocalize. The high PPcoloc demonstrated robust evidence for colocalization between these traits at rs429358.

Finally, a cluster between T2D and T2D adjusted for BMI was identified in the main analysis but was not replicated in the secondary analysis (Table 1). Instead, a cluster between triglycerides and hip circumference adjusted for BMI was identified, driven by an independent variant, rs5117 (R2 with rs8108269 < 0.001) (Table 1). This discrepancy is likely to be a result of the number of variants present in the main analysis.

Conditional Analysis at the GIPR Locus

Our univariate two-sample MR results showed that E354 was associated with a total of 20 traits at a nominal significance threshold (Fig. 1). Independent signal selection showed that E354 or proxy variants in high LD (R2 > 0.9) with E354 were identified as independent signals for fasting GIP, 2-h glucose, total cholesterol levels, BMI, and X-12283 levels. A total of 24 variants were independently associated with CAD on chromosome 19, four of which were in the 1-Mb regions either side of E354 at the GIPR locus (Table 2). Conditioning the association between E354 and CAD on the residual LD between E354 and rs7412, the variant estimated to drive the cluster with CAD, resulted in a slight attenuation of this association but remained significant (OR per copy of E354 after adjustment 1.03; 95% CI 1.02, 1.04; P = 0.003). Of the independent variants identified, rs1964272, an intronic variant in small nuclear ribonucleoprotein D2 polypeptide (SNRPD2), was estimated to be in the strongest LD with E354 (R2 = 0.27) (Fig. 3 and Supplementary Fig. 9). The association between E354 and CAD risk was attenuated when conditioned on rs1964272 (OR per copy of E354 after adjustment 1.01; 95% CI 0.99, 1.03; P = 0.06) (Table 3). In line with this, the association between rs1964272 and CAD risk was attenuated but remained significant with conditioning on E354 (β per copy of rs1964272 after adjustment 0.02; 95% CI 0.01, 0.03; P = 7 × 10−4) (Supplementary Table 6). In addition, the association between E354 and small vessel stroke was also attenuated when conditioned on rs1964272 (Table 3). None of the other loci previously shown to be associated with fasting GIP levels were found to be associated with CAD (Supplementary Table 7). Interestingly, rs1964272 was also associated with levels of QPCTL and SCGB3A1, indicating confounding by LD for the proteomics data as well (Supplementary Fig. 10). Conditioning the association between E354 and QPCTL levels on rs1964272 attenuated the association to nonsignificance (β QPCTL per copy of E354 after adjustment 0.01; 95% CI −0.02, 0.04; P = 0.48) (Table 3).

Table 2

Independent CAD variants identified using approximate conditional analysis

Variant*	Chr: pos	Closest gene	EA	EAF	Marginal β (SE)†	Marginal P†	Conditional β (SE)‡	Conditional P‡	N	R2 with rs1800437
rs429358	19: 45411941	APOE	T	0.85	−0.09 (0.008)	2.86 × 10−27	−0.08 (0.008)	5.87 × 10−23	286,423	0.001
rs7412	19: 45412079	APOE	T	0.08	−0.14 (0.011)	1.66 × 10−35	−0.12 (0.011)	1.58 × 10−28	275,803	0.004
rs11673093	19: 45742094	EXOC3L2	A	0.26	0.04 (0.007)	4.11 × 10−11	0.04 (0.007)	3.09 × 10−10	300,789	0
rs1964272	19: 46190268	SNRPD2	A	0.48	−0.03 (0.006)	9.65 × 10−9	−0.03 (0.006)	1.87 × 10−7	299,519	0.27

Chr, chromosome; EA, effect allele; EAF, effect allele frequency; pos, position; N, number of participants; R2, LD estimate.

The independent CAD variants in the 1-Mb region on either side of E354 are shown.

Log ORs from the original GWAS summary statistics.

Log ORs from the joint model fitted by GCTA.

Figure 3

Regional association plots depicting CAD lead variants in the GIPR region. A: The independent CAD lead variants in the GIPR region are labeled, and their respective associations with CAD are shown before conditional analysis. The region around rs1800437 (E354) is expanded in the red insert to show the LD and proximity of rs1964272 to rs1800437. B: The associations of variants in the GIPR region after conditioning on rs1964272. The region around rs1800437 (E354) is expanded in the red insert to show the attenuation of the E354 signal when conditioned on rs1964272.

Table 3

Conditioning each of the traits associated with E354 at nominal significance from the 2SMR analysis on independent SNPs for each trait

	2SMR result		Conditional result		Independent variant
Trait	β (SE)	P	β (SE)	P ‡	Conditioned on*	LD with rs1800437†
T2D	−0.03 (0.007)	7 × 10−5	−0.03 (0.008)	4 × 10−4	rs3810291	0.001
T2DadjBMI	−0.07 (0.009)	2 × 10−14	−0.02 (0.009)	0.04	rs2238689	0.363
CAD	0.03 (0.007)	2 × 10−6	0.01 (0.007)	0.06	rs1964272	0.269
SVS	−0.08 (0.029)	0.009	−0.04 (0.029)	0.12	rs1964272	0.269
Nonfasted plasma glucose	0.02 (0.003)	3 × 10−8	0.01 (0.003)	0.05	rs1964272	0.269
HbA1c	−0.01 (0.003)	1 × 10−7	−0.0003 (0.003)	0.92	rs9676912	0.356
ApoA1	0.01 (0.003)	3 × 10−6	0.002 (0.003)	0.37	rs2238689	0.363
HDL	0.01 (0.003)	5 × 10−7	0.003 (0.003)	0.25	rs2238689	0.363
ApoB	0.02 (0.002)	5 × 10−13	0.01 (0.002)	2 × 10−5	rs7412	0.004
LDL	0.02 (0.003)	2 × 10−16	0.016 (0.003)	1 × 10−8	rs7412	0.004
Triglycerides	−0.01 (0.003)	2 × 10−5	−0.01 (0.003)	5 × 10−5	rs4803936	0.001
CRP	−0.01 (0.002)	0.02	−0.004 (0.002)	0.07	rs7412	0.004
Albumin	−0.01 (0.003)	6 × 10−6	−0.01 (0.003)	0.001	rs35114617	0.061
Creatinine	−0.02 (0.002)	1 × 10−11	−0.02 (0.002)	3 × 10−11	rs7412	0.004
QPCTL	−0.07 (0.016)	9 × 10−6	0.01 (0.016)	0.48	rs1964272	0.269
Secretoglobin family 3A member 1	−0.08 (0.017)	6 × 10−7	−0.04 (0.017)	0.01	rs61703905	0.1

Estimates of 2-h glucose, total cholesterol, and BMI were not included in this table, as the independent signal selection showed that E354 was one of the independent variants. adj, adjusted for; CRP, C-reactive protein; QPCTL, glutaminyl-peptide cyclotransferase like; SVS, small vessel stroke.

The independent variant showing the greatest attenuation of the E354 association estimate with the respective trait.

LD estimates are in R2 and are quoted from five European populations in the LDlink database v4.1.0.

A nominal significance threshold of P ≤ 0.05 was used to ascertain significance for the conditional results.

Conditioning the association of E354 with LDL, apoB, and triglycerides on independent variants for each trait showed that these remained statistically significant despite being attenuated (Table 3), suggesting that E354 may have independent effects on lipid metabolism.

Discussion

In this study, we applied Bayesian multitrait colocalization and conditional analysis to gain greater understanding of the underlying genetic architecture of CAD and its relation to fasting GIP levels at the GIPR locus. Multitrait colocalization robustly identified a cluster of CAD and lipid traits at APOE that was independent from a cluster of fasting and 2-h GIP, glycemic and anthropometric traits driven by E354. Further, conditional analysis robustly attenuated E354’s association with CAD, small vessel stroke, and QPCTL levels with adjustment for rs1964272 in SNRPD2, an established CAD risk locus (53). Together these results show that association signals for CAD at GIPR are not mediated by an independent effect of GIPR variants on CAD risk but are instead the result of LD confounding between E354 and rs1964272.

Taken together, these findings highlight the specificity of E354’s effects on fasting GIP levels and robustly demonstrate that higher E354-mediated fasting GIP levels are not associated with CVD risk. These results contradict recent genetic evidence linking higher fasting GIP levels with increased CVD risk (21,27), which led to concerns that chronic pharmacological GIPR agonism could have detrimental effects on cardiovascular health (27) and represent safety concerns for pharmacological agonism of this pathway (54). We therefore provide evidence that the inclusion of GIPR agonism in dual GIPR/GLP1R agonists could potentiate the protective effect of GLP-1 agonists on diabetes without undue CVD risk, an aspect not yet assessed in clinical trials. Many studies have shown that GLP1R agonism achieved through chronic pharmacologic therapy, or genetic gain of function, is associated with improved cardiovascular outcomes (28–32). Hence, the available evidence suggests that dual agonism of these receptors may exploit the metabolically favorable combined pharmacology of these incretins without undue CVD risk. However, this proposition requires formal assessment in clinical trials such as the recently initiated cardiovascular outcomes trial SURPASS-CVOT (cardiovascular outcomes trial) of the GIP/GLP1R dual agonist tirzepatide (clinical trial reg. no. NCT04255433, ClinicalTrials.gov).

This study has potential limitations. Firstly, our analysis focuses on a single locus associated with both fasting GIP levels and CAD. This assumes that the GIPR locus is a suitable proxy for fasting GIP levels within which to partition the associations of these two complex traits. Considering that the association at this locus with 2-h glucose is statistically robust and in line with the established function of GIP, this is a reasonable assumption. In addition, no other locus has been reported to be associated with both fasting GIP and CAD, and examining the association of other variants associated with fasting GIP levels (36) in genes other than GIPR showed no association of any of these variants with CAD. However, this does not preclude the existence of other variants that have not yet been associated with GIP levels and may contribute to CVD risk. Patients with T2D are the target of GIPR/GLP1R agonist treatment. We investigate the genetic association of E354 on CAD using the largest publicly available genome-wide summary statistics (53). Therefore, analyses stratified by T2D status are not possible, since such results were not generated and are, hence, not available. Indeed, pursuing this in individual studies would vastly lower sample sizes and therefore be underpowered to detect whether associations with CAD differ significantly by T2D status. Specifically, to affect our results and conclusions about the E354-CAD association being the result of confounding by LD, the genetic architecture at GIPR would have to differ between European-descent individuals with and without prevalent T2D, such that the residual confounding by LD differs by T2D status. As LD is generally preserved between individuals from the same ethnic group, this is a very unlikely scenario.