eQTL mapping in fetal-like pancreatic progenitor cells reveals early developmental insights into diabetes risk

Abstract

The impact of genetic regulatory variation active in early pancreatic development on adult pancreatic disease and traits is not well understood. Here, we generate a panel of 107 fetal-like iPSC-derived pancreatic progenitor cells (iPSC-PPCs) from whole genome-sequenced individuals and identify 4065 genes and 4016 isoforms whose expression and/or alternative splicing are affected by regulatory variation. We integrate eQTLs identified in adult islets and whole pancreas samples, which reveal 1805 eQTL associations that are unique to the fetal-like iPSC-PPCs and 1043 eQTLs that exhibit regulatory plasticity across the fetal-like and adult pancreas tissues. Colocalization with GWAS risk loci for pancreatic diseases and traits show that some putative causal regulatory variants are active only in the fetal-like iPSC-PPCs and likely influence disease by modulating expression of disease-associated genes in early development, while others with regulatory plasticity likely exert their effects in both the fetal and adult pancreas by modulating expression of different disease genes in the two developmental stages.

Fig. 1. Discovery and characterization of eQTLs in iPSC-PPC.

a Study overview created using PowerPoint. b Density plots showing the distribution of PDX1⁺ cells (%; regardless of NKX6-1 status; light green) and PDX1⁺/NKX6-1⁺ cells (%; dark green) in the 107 iPSC-PPC samples. c Bar plot showing the number of eGenes with primary and conditional e_gQTLs. d Bar plot showing the number of eIsoforms with primary and conditional e_iQTLs. e Enrichment (odds ratio, X-axis) of eQTLs in genomic regions (Y-axis) using two-sided Fisher’s Exact Tests comparing the proportion of variants with causal posterior probability (PP) $\ge$5% in the genomic regions between e_gQTLs (blue; n = 8763) and e_iQTLs (yellow; n = 8919). f Line plot showing Pearson correlations of TF binding score and eQTL effect size at different thresholds of causal PP for e_gQTLs (blue) and e_iQTLs (yellow) (Supplementary Data 8). Closed points indicate significant correlations (nominal p < 0.05) while open points indicate non-significant correlations (nominal (p > 0.05). g Bar plot showing the number of genes with only e_gQTLs (blue; n = 3057), only e_iQTLs (green; n = 1554), or both. Orange represents genes whose e_gQTLs colocalized with all their corresponding e_iQTLs (PP.H4 $\ge$ 80%; n = 333). Red represents genes whose e_gQTLs did not colocalize with any of their corresponding e_iQTLs (PP.H3 $\ge$ 80%; n = 38), and pink represents genes with both shared and distinct e_gQTLs and e_iQTLs (i.e., an eGene with two eIsoforms may colocalize with one eIsoform but not the other) (n = 39). Gray represents genes whose eQTL signals were not sufficiently powered to test for colocalization (PP.H4 < 80% and PP.H3 < 80%; n = 598).

Fig. 2. Comparison of the genetic architecture underlying gene expression between fetal-like and adult islets.

a Venn diagram showing the overlap of eGenes between fetal-like iPSC-PPC and adult islets. b Stacked bar plot showing the total number of eGenes detected in adult islets (blue; n = 4211 total) that were expressed in iPSC-PPC (light blue; n = 3301). Likewise, we show the total number of fetal-like iPSC-PPC eGenes (green; n = 4065 total) that were expressed in adult islets (light green; n = 3605). These results show that the majority of eGenes were expressed in both tissues, however, a large fraction was influenced by genetic variation in only one of the two tissues. Therefore, the small overlap of eGenes may be due to differences in the genetic regulatory landscape. c Pie chart showing the proportion of shared eGenes with distinct genetic loci (PP.H3 $\ge$ 80%, purple) or shared genetic loci (PP.H4 $\ge$ 80%, orange). These results show that 12% of the shared eGenes were associated with distinct regulatory variants between fetal-like and adult pancreatic stages. d Example of a shared eGene (SNX29) whose expression was associated with distinct e_gQTL signals (PP.H3 = 90.4%) in fetal-like iPSC-PPC (green, top panel) and adult islets (blue, bottom panel). The X-axis represents variant positions while the Y-axis shows the −log10(eQTL p-value) for the associations between the genotype of the tested variants and gene expression. For plotting purposes, we assign a single p-value for gene-level significance after Bonferroni-correction (0.05/number of independent variants tested in fetal-like iPSC-PPC; horizontal line). Red vertical lines show the positions of the lead variants in fetal-like iPSC-PPC and adult islets (chr16:12656135:C > G and chr16:12136526:A > G, respectively).

Fig. 3. eQTL sharing between iPSC-PPC, adult islets, and adult whole pancreas.

a Bar plot showing the number of tissue-unique e_gQTLs identified in fetal-like iPSC-PPC, adult islets and adult whole pancreas. b Bar plot showing the number of e_gQTL modules for each annotation. c Top panels: Enrichment (odds ratio) of iPSC-PPC singleton and combinatorial e_gQTLs in hESC-derived PPC chromatin states. Bottom panels: Enrichment (odds ratio) of iPSC-PPC singleton and combinatorial e_gQTLs in adult islet chromatin states. Enrichment was calculated using a two-sided Fisher’s Exact Test comparing the proportion of candidate causal variants overlapping the chromatin states versus a background of randomly selected 20,000 variants at various PP thresholds. P-values were Benjamini-Hochberg-corrected and considered significant if the corrected p-values < 0.05 (indicated by asterisk, Supplementary Data 12). d CDC37L1-DT locus showing an iPSC-PPC-unique singleton e_gQTL overlapping an adult islet active promoter region. Lower panel shows the positions of active promoters in the adult islets. e, f The chr3:148903264-148983264 locus (gray rectangle) is an example of an “iPSC-PPC-unique” module (module ID: GE_3_1) associated with CP and HPS3 expression. g, h The chr15:57746360-57916360 locus (gray rectangle) is as an example of an “adult islet-unique” module (module ID: GE_15_13) associated with GCOM1, MYZAP, and POLR2M expression. GCOM1 was not expressed in adult whole pancreas and therefore, was not tested for e_gQTL association. i, j The chr5:146546063-146746063 locus (gray rectangle) is an “adult whole pancreas-unique” e_gQTL module (module ID: GE_5_32) associated with STK32A and STK32A-AS1 expression. STK32A-AS1 was not expressed in iPSC-PPC and therefore, was not tested for e_gQTL association. Panels e, g, i display the e_gQTL modules as networks in which the e_gQTL associations (nodes) are connected by edges due to colocalization (PP.H4 $\ge$ 80%). For panels d, f, h, and j, the X-axis represents variant positions while the Y-axis shows the −log10(eQTL p-value) for the associations between the genotype of the tested variants and gene expression. For plotting purposes, we assigned a single p-value for gene-level significance after Bonferroni-correction (0.05/the number of independent variants tested in fetal-like iPSC-PPC; horizontal line). Red vertical lines indicate the positions of the lead candidate causal variants underlying the colocalization based on maximum PP.

Fig. 4. Regulatory plasticity of e_gQTL loci.

a Number of e_gQTL modules shared between iPSC-PPC and at least one adult pancreas tissue categorized by eGene overlap with adult. “Zero” indicates that the module did not contain an e_gQTL in the respective adult tissue. “Same” indicates that the module had e_gQTLs for only the same eGenes in iPSC-PPC and the adult tissue. “Partial” indicates that the module had e_gQTLs for partially overlapping eGenes between iPSC-PPC and the adult tissue. “Different” indicates that the module had e_gQTLs for only different eGenes between iPSC-PPC and the adult tissue. b–d Examples of e_gQTL loci demonstrating regulatory plasticity of genetic variation across fetal-like and adult pancreatic stages. Panel b shows a locus strongly associated with AC119427.1 expression in fetal-like iPSC-PPC and TNNI1 expression in adult islet and whole pancreas. Panel c shows a locus associated with MPND expression in only fetal-like iPSC-PPC but STAP2 expression in both the adult pancreatic tissues. Panel d shows a locus associated with partially overlapping eGenes between the two pancreatic stages (UROS in all three pancreatic tissues and BCCIP in only adult islets). The X-axis represents variant positions while the Y-axis shows the −log10(eQTL p-value) for the associations between the genotype of the tested variants and gene expression. For plotting purposes, we assigned a single p-value for gene-level significance after Bonferroni-correction (0.05/the number of independent variants tested in fetal-like iPSC-PPC; horizontal line). Red vertical lines indicate the positions of the lead candidate causal variants underlying the colocalization based on maximum PP.

Fig. 5. Summary of pancreatic GWAS associations.

a Bar plot showing the number of eQTL loci that colocalized with GWAS variants (PP.H4 $\ge$ 80%) as a singleton or module. b Pie chart showing the number of singleton-colocalized GWAS loci (n = 183) color-coded by the number of candidate causal variants identified in their 99% credible sets. c Pie chart showing the number of module-colocalized GWAS loci (n = 129) color-coded by the number of candidate causal variants identified in their 99% credible sets.

Fig. 6. Pancreatic GWAS associations with fetal-specific and adult-shared gene Expression.

a The TPD52 locus is associated with fasting glucose levels and colocalized with a fetal-like iPSC-PPC-unique singleton e_gQTL with the predicted causal variant identified as rs12549167 (chr8:81078464:C > T, PP = 33.9%). b The CDC37L1-DT locus is associated with fasting glucose and type 1 diabetes and colocalized with an iPSC-PPC-unique singleton e_gQTL with the predicted causal variant identified as rs10758593 (chr9:4292083:G > A, PP = 79.2%). c Cholesterol and LDL direct GWAS loci colocalize with a fetal-adult e_gQTL module where the variants are strongly associated with ADSL expression in iPSC-PPC and ST13 expression in the adult whole pancreas (also weakly associated with ST13 expression in the adult islets). The predicted causal variant was identified as rs138349 (chr22:41249522:A > G, PP = 21.9%). For plotting purposes, we assigned a single p-value for gene-level significance based on Bonferroni-correction (0.05 divided by the number of independent variants tested in fetal-like iPSC-PPC; horizontal line). We note that this p-value does not reflect the thresholds used to define a significant eQTL in the original adult studies^,. Therefore, while the ST13 eQTL in adult islets in panel c is below the horizontal line, it had an FDR < 1% in the original study. In each panel, the X-axis represents variant positions while the Y-axis either shows the −log10(eQTL p-value) for the associations between the genotype of the tested variants and gene expression or the −log10(GWAS p-value) for the associations between the tested variants and the GWAS trait. For GWAS significance, we used −log10(5 × 10⁻⁸). Red vertical lines indicate the positions of the lead candidate causal variants underlying the colocalization based on maximum PP. For loci that colocalized with multiple GWAS traits, we used the credible set that yielded the smallest number of variants to plot the “PP” fine-mapping panel.

Fig. 7. Pancreatic GWAS associations with fetal-specific alternative splicing.

a T1D-risk locus colocalized with a fetal-like iPSC-PPC-unique singleton e_ASQTL for MEG3 with the predicted causal variant rs56994090 (chr14:101306447:T > C, PP = 100%). b GWAS locus associated with HbA1c colocalized with an iPSC-PPC-unique singleton e_ASQTL for CDH3 with the predicted causal variant rs72785165 (chr16:68755635:T > A, PP = 6.8%). c GWAS locus associated with T2D-risk and BMI colocalized with an iPSC-PPC-unique e_ASQTL module (AS_13_2) for differential usage of three HMGB1 isoforms with a predicted causal variant rs3742305 (chr13:31036642:C > G, PP = 49.3%). In each panel, the X-axis represents variant positions while the Y-axis either shows the −log10(eQTL p-value) for the associations between the genotype of the tested variants and gene expression or the −log10(GWAS p-value) for the associations between the tested variants and the GWAS trait. For GWAS significance, we used −log10(5 × 10⁻⁸). For eQTL significance, we used a single p-value for gene-level significance after Bonferroni-correction (0.05/the number of independent variants tested in fetal-like iPSC-PPC; horizontal line). Red vertical lines indicate the positions of the lead candidate causal variants underlying the colocalization based on maximum PP. For loci that colocalized with multiple GWAS traits, we used the credible set that yielded the smallest number of variants to plot the “PP” fine-mapping panel.