Adipose tissue eQTL meta-analysis reveals the contribution of allelic heterogeneity to gene expression regulation and cardiometabolic traits

Abstract

Complete characterization of the genetic effects on gene expression is needed to elucidate tissue biology and the etiology of complex traits. Here, we analyzed 2,344 subcutaneous adipose tissue samples and identified 34K conditionally distinct expression quantitative trait locus (eQTL) signals in 18K genes. Over half of eQTL genes exhibited at least two eQTL signals. Compared to primary signals, non-primary signals had lower effect sizes, lower minor allele frequencies, and less promoter enrichment; they corresponded to genes with higher heritability and higher tolerance for loss of function. Colocalization of eQTL with conditionally distinct genome-wide association study signals for 28 cardiometabolic traits identified 3,605 eQTL signals for 1,861 genes. Inclusion of non-primary eQTL signals increased colocalized signals by 46%. Among 30 genes with ≥2 pairs of colocalized signals, 21 showed a mediating gene dosage effect on the trait. Thus, expanded eQTL identification reveals more mechanisms underlying complex traits and improves understanding of the complexity of gene expression regulation.

Figure 1.. Conditionally distinct signals in adipose eQTL studies.

A. Number of genes with 1 to 10 eQTL signals (P ≤ 1e-6) identified in each study and the meta-analyses. ‘≥1 signal’ column indicates the number of genes with at least one significant eQTL signal, ‘≥2 signals’ indicates the number of genes with two or more eQTL signals, and the percentage of genes with an eQTL that have two or more eQTL signals is in parentheses. B. The numbers of genes identified with an eQTL in each study are represented by filled circles, and the numbers of eQTL signals are represented by asterisks. Studies are shown by color: blue, METSIM (N); purple, METSIM (S); green, FUSION; orange, TwinsUK; pink, GTEx all populations; red, GTEx EUR; gray, meta-analysis with GTEx all populations; and black, meta-analysis with GTEx EUR. C. The number of genes with 1 through 10 eQTL signals detected in each study.

Figure 2.. GLYCTK eQTL signals identified in each study and the meta-analysis.

LocusZoom plots of the marginal GLYCTK eQTL for the meta-analysis and each individual study. The x-axes show position on chromosome 3 and y-axes show eQTL −log₁₀ P-value. The lead variant of the 1^st signal (chr3:52,273,421, rs610060) in the meta-analysis is represented by a red diamond in all plots, and the lead variant of the second signal (chr3:52,276,901, rs11711914) in the meta-analysis is represented by a blue diamond in all plots. The red circles represent variants in stronger LD with the lead variant of the 1^st signal while the blue triangles represent variants in stronger LD with the lead variant of the second signal. Shading indicates LD r² as shown in the legend. Although each study has both signals colored, only one signal was significant in the conditional eQTL analysis for each of the individual studies (P ≤ 1e-6).

Figure 3.. Characteristics of eQTL variants and genes according to the number of significant eQTL signals.

Violin plots with inset boxplots of the (A) absolute value of the effect sizes of lead variants, (B) MAF, and (C) distance of the lead variants to the gene TSS for the indicated signals in order of discovery. Only the 661 genes with 5 or more signals were included. For the boxplots, the center line represents the median value, the box limits represent the upper and lower quartiles, whiskers represent the 1.5x interquartile range, and the black circles represent outliers. The black lines connect the median values of each signal group. In C, 163 points with a distance to TSS greater than 600 were excluded. See Figure S6 for genes with one to four eQTL signals. (D) Proportion of genes in TwinsUK with the specified number of eQTL signals separated by gene expression quartiles. Quartile 1 indicates the genes with the lowest expression. The darkest blue are the genes without an eQTL signal and the lightest blue are genes with five or more eQTL signals. (E) Violin plots with inset boxplots of the heritability of genes with the specified number of eQTL signals in TwinsUK. For the boxplots, the center line represents the median value, the box limits represent the upper and lower quartiles, whiskers represent the 1.5x interquartile range, and the black circles represent outliers. The black lines connect the median values of each signal group. (F) Proportion of genes for each signal number with a pLI score ≥ 0.9 out of the total number of genes that have pLI scores available for that signal number.

Figure 4.. Sex-stratified WHRadjBMI GWAS and ADORA1 eQTL signal plots.

(A) LocusZoom plots for WHRadjBMI female GWAS signal and (B) ADORA1 female eQTL signal. (C) LocusZoom plots for WHRadjBMI male GWAS signal and (D) ADORA1 male eQTL signal. All plots are colored by LD with the female GWAS lead variant represented by a purple diamond.

Figure 5.. Colocalization of two or more GWAS signals with two or more eQTL signals at ZNRF3 and PDE3A.

A. LocusZoom plots of WHRadjBMI GWAS summary statistics (Pulit et al 2019) (top) and marginal ZNRF3 eQTL data for the meta-analysis (bottom). Both plots show two signals colored by the GWAS lead variants (red diamond, 1^st signal chr22:29,449,477, rs2294239; blue diamond, 2^nd signal chr22:29,338,235, rs5762906). The red circles and blue triangles indicate genetic variants in stronger LD with the 1^st or 2^nd signal, respectively and are shaded based on LD. Signal 1 in the GWAS is colocalized with signal 1 of the eQTL dataset (LD r² = 0.90; coloc PP4 = 0.99) and signal 2 for both datasets are also colocalized (LD r² = 1.00; coloc PP4 = 0.98). B. Effect sizes of the WHRadjBMI GWAS signals (y-axis) versus the effect sizes of the ZNRF3 eQTL signals (x-axis) from MRLocus. Each point represents a colocalized eQTL signal with standard error bars. The solid blue line represents the slope of the effect of the gene on the trait, and dotted blue lines represent the confidence interval. The slope estimates a gene-to-trait effect of −0.19, meaning that increasing adipose ZNRF3 expression level by one population standard deviation should reduce WHRadjBMI by 19% of its population standard deviation. C. Scatter plot of inverse normalized waist-to-hip ratio (x-axis) and ZNRF3 gene expression (y-axis) in METSIM (S) (n = 420). Each point represents an individual sample, the blue line represents the linear regression slope and the 95% confidence interval is shown in gray. The correlation value and association P-value are shown. D. LocusZoom plot of the HDL-C GWAS summary statistics (Graham et al 2021) (top) and marginal PDE3A eQTL data for the meta-analysis at (bottom). Both plots show four signals colored by the GWAS lead variants (red diamond, 1^st signal chr12:20,470,221, rs11045172; blue diamond, 2^nd signal chr12:20,470,009, rs2044315; yellow diamond, 3^rd signal chr12:20,579,083, rs11045237; green diamond, 4^th signal chr12:20,591,332, rs7134150). The red circles, blue triangles, yellow squares, and green inverted triangles indicate genetic variants in stronger LD with the 1^st, 2^nd, 3^rd, or 4^th signal, respectively and are shaded based on LD. Signal 1 in the GWAS is colocalized with signal 1 of the eQTL dataset (LD r² = 1.00; coloc PP4 = 1.00), signal 2 for the GWAS is colocalized with signal 4 of the eQTL dataset (LD r² = 0.93; coloc PP4 = 1.00), signal 3 for the GWAS and signal 2 for the eQTL dataset are colocalized (LD r² = 0.42; coloc PP4 = 1.00), and signal 4 for the GWAS and signal 3 for the eQTL dataset are colocalized (LD r² = 0.94; coloc PP4 = 0.99). E. Effect sizes of the HDL-C GWAS signals (y-axis) versus the effect sizes of the PDE3A eQTL signals (x-axis) from MRLocus. Each point represents a colocalized eQTL signal with standard error bars. The solid blue line represents the slope of the effect of the gene on the trait, and dotted blue lines represent the confidence interval. The slope estimates a gene-to-trait effect of −0.14, meaning that increasing adipose PDE3A expression level by one population standard deviation should reduce HDL-C by 14% of its population standard deviation. F. Scatter plot of inverse normalized HDL-C (x-axis) and PDE3A gene expression (y-axis) in METSIM (S) (n = 420). Each point represents an individual sample, the blue line represents the linear regression slope and the 95% confidence interval is shown in gray. The correlation value and association P-value are shown.

Figure 6.. Regulatory annotation enrichment of eQTL signals and validation of allelic effects on transcriptional activity at SEMA3C.

A. eQTL signals enriched in Roadmap Epigenomics chromatin states in adipose tissue compared to genes without an eQTL separated by signal number. Dark red represents promoters and gold represents enhancers. The bars represent the upper and lower 95% confidence intervals. The asterisk represents significant Bonferroni-adjusted enrichment values that do not overlap an odds ratio (OR) of 1 (black dashed line). B. LocusZoom plots of the WHRadjBMI GWAS summary statistics (Pulit et al 2019)(top) and the SEMA3C meta-analysis eQTL data conditioned on all but signal 1 (bottom). Both plots show the same lead variant represented by a purple diamond (chr7:80,570,871; rs917191). Other variants are colored based on the LD r² 1000G EUR with the lead variant. Signal 1 in the GWAS dataset is colocalized with signal 1 of the eQTL dataset (LD r² = 1.0; coloc PP4 = 1.0). C. UCSC genome browser tracks showing regulatory annotations that overlap SEMA3C eQTL variants. In the SEMA3C SNPs track, the lead variant is shown in purple and proxy variants (LD r² ≥ 0.8) are shown in black. The chromHMM tracks are from Epigenomic Roadmap for mesenchymal stem cell-derived adipocytes, adipose nuclei, skeletal muscle, liver, and brain hippocampus; red represents a promoter-like signature, yellow represents an enhancer-like signature, green represents a signature for elongating RNA polymerase, and gray represents low to no signal. The blue signal tracks represent ATAC-seq accessible chromatin in SGBS cells at differentiation day 0, day 4, and day 14. The METSIM adipose peaks are ATAC-seq peaks detected in at least 3 adipose tissue samples. SEMA3C gene annotations are from UCSC genes. The bottom figure shows the browser tracks zoomed in to the region around rs917191. D. Relative transcriptional activity of rs917191-G and rs917191-C in hWAT adipocytes from dual-luciferase reporter assays. Values indicate transcriptional activity relative to an empty vector (EV), points represent independent clones with standard error bars, and P-values from Student’s unpaired t-tests compare activity between alleles. E. Scatter plot of inverse normalized waist-to-hip ratio (y-axis) and SEMA3C gene expression (x-axis) in METSIM (S) (n = 420). Each point represents an individual sample, the blue line represents the linear regression slope and the 95% confidence interval is shown in gray. The correlation value and association P-value are shown.