GWAS of random glucose in 476,326 individuals provide insights into diabetes pathophysiology, complications and treatment stratification

Abstract

Conventional measurements of fasting and postprandial blood glucose levels investigated in genome-wide association studies (GWAS) cannot capture the effects of DNA variability on 'around the clock' glucoregulatory processes. Here we show that GWAS meta-analysis of glucose measurements under nonstandardized conditions (random glucose (RG)) in 476,326 individuals of diverse ancestries and without diabetes enables locus discovery and innovative pathophysiological observations. We discovered 120 RG loci represented by 150 distinct signals, including 13 with sex-dimorphic effects, two cross-ancestry and seven rare frequency signals. Of these, 44 loci are new for glycemic traits. Regulatory, glycosylation and metagenomic annotations highlight ileum and colon tissues, indicating an underappreciated role of the gastrointestinal tract in controlling blood glucose. Functional follow-up and molecular dynamics simulations of lower frequency coding variants in glucagon-like peptide-1 receptor (GLP1R), a type 2 diabetes treatment target, reveal that optimal selection of GLP-1R agonist therapy will benefit from tailored genetic stratification. We also provide evidence from Mendelian randomization that lung function is modulated by blood glucose and that pulmonary dysfunction is a diabetes complication. Our investigation yields new insights into the biology of glucose regulation, diabetes complications and pathways for treatment stratification.

Fig. 1. Summary of all RG loci identified in this study.

a, Circular Manhattan plot summarizing findings from this study. In the outermost layer, gene names of the 133 distinct RG signals are labeled with different colors indicating the following three clusters defined in cluster analysis: 1a/1b, metabolic syndrome; 2a/2b, insulin release versus insulin action (with additional effects on inflammatory bowel disease for cluster 2a) and 3, defects of insulin secretion. Asterisks annotate RG signals that are new for glycemic traits. Track 1 shows RG Manhattan plot reporting −log₁₀(P value) for RG GWAS meta-analysis. Signals reaching genome-wide significance (P < 5.0 × 10⁻⁸) are colored in red. Crosses annotate loci that show evidence of sex heterogeneity (P_{sex-dimorphic} < 5.0 × 10⁻⁸ and P_{sex-heterogeneity} < 0.05); blue crosses for larger effects in men, green crosses for larger effects in women. Track 2 shows the effects of the 133 independent RG signals on four GIP/GLP-1-related traits GWAS. The colors of the dotted lines indicate four GIP/GLP-1-related traits: gray dot, signals reaching P < 0.010 for a GIP/GLP-1-related trait; red dot, lead SNP has a significant effect on GIP/GLP-1-related trait (Bonferroni corrected P < 1.0 × 10⁻⁴). Track 3 shows the effects (−log₁₀(P value)) of the 133 independent RG signals on 113 glycan PheWAS. Track 4 shows the effects (−log₁₀(P value)) of the 133 independent RG signals on 210 gut-microbiome PheWAS. Track 5 shows MetaXcan results for ten selected tissues for RG GWAS meta-analysis; signals colocalizing with genes (Bonferroni corrected P < 9.0 × 10⁻⁷) are plotted for each tissue. All P values were calculated from the two-sided z statistics computed by dividing the estimated coefficients by the estimated standard error, without adjustment. b, Credible set analysis of RG associations in the European ancestry meta-analysis. Variants from each of the RG signal credible sets are grouped based on their posterior probability (the percentiles labeled on the sides of the bar). SNP variants with posterior probability >80%, along with their locus names, are provided. All variants from the credible set of lead signals are highlighted in bold.

Fig. 2. Functional and structural analysis of coding GLP1R variants.

a, Minor allele frequency-weighted linear regression was used to test if mini-G_s response to GLP-1 stimulation substantially predicted point estimates of GLP1R variant effect on RG levels (AST20 β_RG as estimated in the UKBB study, n_max = 401,810). Mini-G_s response to GLP-1 stimulation was corrected for variant surface expression (n_max = 22, exact n for each variant is provided in Supplementary Table 11). Error bars extend one standard error above and below the point estimate. Size of the dots is proportional to the weight applied in the regression model. The regression results (coefficient of determination R² = 0.74, F(1, 15) = 47.5, P = 5.1 × 10⁻⁶) suggest that mini-G_s coupling in response to GLP-1 stimulation predicts the effect of these coding variants on RG levels (AST20 β_RG = −0.030; 95% confidence interval (CI) = −0.039 to −0.020; P = 5.1 × 10⁻⁶). The gray shaded area around the regression line corresponds to the 95% CI of predictions from the model. Variants in red showed no detectable surface expression (NDE) and are not included in regression analysis. b, Mean GLP1R variant mini-G_s coupling and receptor endocytosis, with surface expression correction, in response to GLP-1, OXM, glucagon (GCG), exendin-4 (Ex4), semaglutide (Sema) and tirzepatide (TZP), n = 6. Positive deviation indicates variant gain-of-function, with statistical significance inferred when the 95% CIs shown do not cross zero. Responses are also compared between pathways by unpaired t test, with an asterisk indicating statistically significant differences. c, Architecture of the complex formed between the agonist-bound GLP-1R and G_s; the likely effect triggered by residues involved in GLP-1R isoforms A316T, G168S and R421W (in magenta) are reported. d, Distributions of the distance between Y242^3.45 side chain and P312^5.42 backbone computed during molecular dynamics simulations of GLP-1R WT and A316T; the cut-off distance for hydrogen bond is shown. e, Difference in the hydrogen bond network between GLP-1R WT and A316T. f, Analysis of water molecules within the TMD of GLP-1R WT and A316T suggests minor changes in the local hydration of position 5.46 (unperturbed structural water molecule). Also, a stabilizing role for the water molecules at the binding site of the G protein (water cluster apha5) cannot be ruled out. g, Distributions of the distance between position 168^1.63 and Y178^2.48 during molecular dynamics simulations of GLP-1R WT and G168S. h, During molecular dynamics simulations, the GLP-1R isoform S168G showed increased flexibility of ICL1 and H8 compared to WT, suggesting a different influence on G-protein intermediate states. i, Contact differences between G_s and GLP-1R WT or W421R; the C terminal of W421R H8 made more interactions with the N terminal segment of the G_s β subunit. j, Mini-G_s and GLP-1R endocytosis responses to 20 nM exendin-4, plotted against surface GLP-1R expression, from 196 missense GLP1R variants transiently transfected in HEK293T cells (n = 5 repeats per assay), with data represented as mean ± s.e.m. after normalization to WT response and log₁₀-transformation. Variants are categorized as ‘LoF1’ when the response 95% CI falls below zero or ‘LoF2’ where the expression-normalized 95% CI falls below zero. k, GLP-1R snake plot created using gpcr.com summarizing the functional impact of missense variants; for residues with >1 variant, classification is applied as LoF2 > LoF1 > tolerated.

Fig. 3. Deterioration of glucose homeostasis progressing into T2D and leading to complications in multiple organs and tissues.

Established (left, in peach) and new (right, in green). a, A human figure illustrating the main causes of hyperglycemia (a combination of lifestyle and genetic factors), and how hyperglycemia affects many organs and tissues. Complications on the left panel are well-established for T2D. Those on the right panel are emerging ones and are supported by our current analyses. Figure created with BioRender.com. b, DEPICT prioritization of 134 tissues from the GTEx Project highlights the ileum and pancreas (shown in red, one-sided empirical P value with FDR < 0.05 determined against randomized phenotypes in a null GWAS).

Fig. 4. Cell type prioritization across 17 tissues identified large intestinal tissue ranked second only to pancreatic cell types.

CELLECT prioritization of 115 cell types from Tabula Muris highlights pancreatic polypeptide (PP) cells (shown in black, one-sided Wilcoxon rank-sum test with significance threshold depicted by a dotted line indicating cell types with a nominal P_S-LDSC < 4.3 × 10⁻⁴).

Fig. 5. Genome-wide genetic correlation between RG and a range of traits and diseases.

The x axis provides the estimated r_g genetic correlation values for traits or diseases (y axis) reaching at least nominal significance (P < 0.05). Correlations reaching P < 0.010 are labeled with the prime symbol, and those P < 2.1 × 10⁻⁴ are labeled with the asterisk symbol. P values were calculated from the two-sided z statistics computed by dividing the estimated r_g by the estimated standard error, without adjustment. Each error bar represents the standard error of the estimate.

Extended Data Fig. 1. RG trait models tested and sensitivity plots showing the correlations between association analyses beta coefficients and Z-scores from RG models in UKBB.

a, The models were labeled according to covariates included and RG cut-offs used. Individuals were included based on two RG cut-offs: <20 mmol/l to account for the effect of extreme RG values (20) and <11.1 mmol/l (11), which is an established threshold for T2D diagnosis. Hence, model 1 – AS20 refers to adjustment for age and sex, using a cut-off of <20 mmol/l, and so forth. b-e, For c, 4,138 individuals were excluded based on HbA1c ≥ 6.5%, in addition to the self-reported or diagnosed T2D cases. Variants with a heterogeneity P-value ≤ 0.05 (beta-coefficient plot) or a Z-score difference between the two models compared >3 (Z-score plots) are annotated. f, An enrichment plot showing the effect of RG signals (AS20 + AST20 model) on T2D. RG and T2D effect sizes are plotted along the y- and x-axes, respectively. Point size is proportional to the statistical significance of the variant for T2D, with red color indicating previously established signals and blue novel signals, respectively. The dashed line represents the line of best fit. Variants with T2D P-value in the lowest decile are labeled. g, An enrichment plot showing the effects of RG signals (AS20 + AST20 model) on HOMA-B and HOMA-IR. The effect sizes on HOMA-B and HOMA-IR are plotted along the y- and x-axes, respectively. Point size is proportional to the significance of the variant either in HOMA-B or HOMA-IR, depending on which trait has the smaller P value. Red color indicates previously established signals and blue indicates novel signals, respectively. Variants with suggestive significance (P < 5.0 × 10⁻⁶) are labeled.

Extended Data Fig. 2. Enrichment plots showing the effect of RG signals (AS20 + AST20 model) on glycemic and respiratory-related phenotypes.

a-i, Look-up of effects was done in previously published genome-wide association studies for HbA1c (a), fasting glucose (b), fasting insulin (c), type 2 diabetes (d), forced expiratory volume in one second (FEV1) (e), forced vital capacity (FVC) (f), FEV1/FVC (g), lung cancer (h) and squamous cell lung cancer (i). RG and other phenotype effect sizes are plotted along the y- and x-axes, respectively. Point size and color are proportional to the significance of the variant in each phenotype, with red indicating higher and blue lower significance, respectively. The dashed line represents the line of best fit. P < 5.0 × 10⁻⁸ was considered statistically significant after adjusting for multiple testing. Two-tailed P-values are reported. Variants with P-values in the lowest decile are labeled.

Extended Data Fig. 3. LocusZoom plots of common variants in UKBB (Europeans) meta-analysis for RG.

a-h, Plots are shown for GCKR (a), TET2 (b), RREB1 (c), NMT1 (d) and WIPI1 (e) loci and low-frequency coding variants at EDEM3 (f), NEUROD1 (g) and GLP1R (h) loci. The x-axis shows the chromosomal position, and the y-axis shows the uncorrected two-sided −log₁₀ P values from the UKBB GWAS conducted using linear mixed-modeling in BOLT-LMM. Horizontal line corresponds to P = 5 × 10⁻⁸ and blue peaks show the recombination rate.

Extended Data Fig. 4. Association analysis of GLPR1 receptor function and random glucose effects of coding variants.

Minor allele frequency-weighted linear regression was used to test if mini-G_s response to GLP-1 stimulation significantly predicted point estimates of GLP1R variant effect on RG levels (AST20 β_RG as estimated in whole-exome sequencing data from the UKBB study). Mini-G_s response to GLP-1 stimulation was corrected for variant surface expression (n_max = 22, exact n for each variant is provided in Supplementary Table 11). Error bars extend one standard error above and below the point estimate. Size of the dots is proportional to the weight applied in the regression model (Methods). The regression results (coefficient of determination R² = 0.56, F(1, 14) = 20.1, P = 5.2 × 10⁻⁴) suggest that mini-G_s coupling in response to GLP-1 stimulation predicts the effect of these coding variants on RG levels (AST20 β_RG = − 0.028; 95% CI = −0.042 to −0.015; P = 5.2 × 10⁻⁴). The gray shaded area around the regression line corresponds to the 95% confidence interval of predictions from the model. Variants in red showed no detectable surface expression (NDE) and are not included in regression analysis.

Extended Data Fig. 5. Epigenetic annotation of the RG GWAS results using GARFIELD.

The analyses were performed using generalized linear modeling in GARFIELD software. We considered enrichment to be statistically significant if the RG GWAS P-value reached P = 1 × 10⁻⁸ and the enrichment analysis P-value was < 2.5 × 10⁻⁵ (Bonferonni corrected for 2,040 annotations).

Extended Data Fig. 6. Cluster analysis of effects (as Z-scores) of the distinct 143 RG signals on 45 relevant phenotypes.

All variant effects were aligned to the RG risk allele. HapMap2 based summary statistics were imputed using SS-Imp v0.5.565 to minimize missingness. Missing summary statistics values were imputed via mean imputation. The heatmap was produced using the Pheatmap package. For visualization, the Z-scores were truncated to the value corresponding to genome-wide significance (Z = 5.45), and 11 phenotypes with the lowest median absolute Z-scores were excluded.

Extended Data Fig. 7. Scatter plots of the standardized allelic effect estimates for selected trait pairs.

In each scatter plot, loci were assigned to the groups defined from the cluster analysis and highlighted by different colors. a, Corrected insulin response (CIR) vs. type 2 diabetes (T2D) (clusters 1a/b related to metabolic syndrome). b, Glycated hemoglobin (HbA1c) vs. inflammatory bowel disease (IBD) (cluster 2a) highlights the effects of loci with a protective role in IBD. c, Plasminogen activator inhibitor-1 (PAI-1) vs. CIR (cluster 3) highlights loci linked to insulin secretion defects.