Identification of an imprinted master trans regulator at the KLF14 locus related to multiple metabolic phenotypes.

Genome-wide association studies have identified many genetic variants associated with complex traits. However, at only a minority of loci have the molecular mechanisms mediating these associations been characterized. In parallel, whereas cis regulatory patterns of gene expression have been extensively explored, the identification of trans regulatory effects in humans has attracted less attention. Here we show that the type 2 diabetes and high-density lipoprotein cholesterol-associated cis-acting expression quantitative trait locus (eQTL) of the maternally expressed transcription factor KLF14 acts as a master trans regulator of adipose gene expression. Expression levels of genes regulated by this trans-eQTL are highly correlated with concurrently measured metabolic traits, and a subset of the trans-regulated genes harbor variants directly associated with metabolic phenotypes. This trans-eQTL network provides a mechanistic understanding of the effect of the KLF14 locus on metabolic disease risk and offers a potential model for other complex traits.

Variants near the maternally-expressed transcription factor KLF14 (Kruppel-like factor 14) are robustly associated with both Type 2 Diabetes (T2D) and HDL-cholesterol levels in large-scale genome-wide association studies (GWAS)1,2. These studies have implicated a group of highly-correlated SNPs including rs4731702 and rs972283 ~14kb upstream of KLF141,2. KLF14 is the regional gene most likely to be mediating these effects since the same SNPs show adipose-specific, maternally-restricted cis-regulatory associations with KLF14 expression levels, a pattern which mirrors the parent-of-origin effects for T2D-susceptibility at this locus3.

Since transcription factors such as KLF14 typically modulate expression of other genes in trans, we tested for association between rs4731702 and expression levels of ~24K probes (16,663 genes) on the Illumina Human HT12 array in subcutaneous adipose tissue biopsies from a cohort of 776 healthy female twins4. The enrichment of rs4731702 trans-associations for low p-values (Figure 1, Supplemental Figure 1) suggests that KLF14 is a master-regulator of gene expression in adipose tissue. The pattern of trans-associations at KLF14 mirrors the GWAS associations (Fig 2, Suppl Fig 2), and conditioning the trans associations on rs4731702 abolishes the signal at all other SNPs. These findings indicate that the same set of SNPs (and presumably the same causal variant) underlies the cis-, trans- and metabolic trait-associations at this locus.

Figure 1

KLF14 is a master regulator of gene expression in adipose tissue

p-value distribution of association between the KLF14 cis-eQTL rs4731702 and expression levels of ~24K probes in adipose tissue.

Figure 2

Regional signal plots of the KLF14 locus

Left: signal plot of a representative trans -regulated gene in the MuTHER samples (N = 776). Middle: signal plot for the Type 2 Diabetes GWA meta-analysis performed by the DIAGRAM consortium1 (stage 1 data only, effective sample size 22,044: rs972283 reached genome-wide significance after further replication). Right: signal plot for the HDL-cholesterol GWA meta-analysis performed by the Lipid Consortium (N = 99,900)2. Signal plots of all 10 genome-wide significant trans-regulated genes are included in Supplementary Figure 2

We focused on the ten genes (TPMT, ARSD, SLC7A10, C8orf82, APH1B, PRMT2, NINJ2, KLF13, GNB1, MYL5) showing genome-wide significant trans (GWST) associations (p < 5×10−8) driven by rs4731702. First, we sought replication of the trans-associations in an independent set of adipose tissue samples (deCODE Genetics; N= 589)5. As previously reported3, the deCODE data revealed a strong maternally-specific cis-association between rs4731702 and KLF14 expression in adipose tissue (p = 1×10−19) (Table1). This cis effect was not detected in the MuTHER data due to apparent problems with the KLF14 probe represented on the Illumina HT12 array used for MuTHER (See methods). Seven of the GWST genes from the MuTHER analysis had a directionally-consistent trans-association with p-value <0.05 in the deCODE replication set (Table 1), and we were able to show parent-of-origin effects for the trans-associations consistent with the maternally-specific cis-effects for KLF14 expression and T2D-risk3. In the deCODE replication data maternally inherited trans-associations were markedly more significant than general analyses and no paternally-inherited trans-associations were seen (Table 1).

Table 1

Genome-wide significant (p< 5 ×10−8) associations of gene expression levels with rs4731702 at 130,083,924 (build 36) on chromosome 7. The effect allele is the Type 2 Diabetes risk allele C, which has a frequency of 55% in the HapMap CEU population.

		MuTHER Data		Decode All		DecodeMaternal		DecodePaternal		Combined MuTHER +Decode Maternal
Gene	Chr	Effect(se)	Pvalue	Effect(se)	Pvalue	Effect(se)	Pvalue	Effect(se)	Pvalue	Zscore	Pvalue	Direction
APH1B	15	0.08(0.013)	6.2E-10	0.11(0.059)	0.08	0.17(0.085)	0.07(0.083)	0.07	0.44	6.2	5.5E-10	+ +
ARSD	X	0.08(0.012)	0.24(0.059)	9.4E-12	2.2E-04	0.51(0.083)	2.6E-08	−0.004(0.083)	0.96	8.7	2.6E-18	+ +
C8orf82	8	0.09(0.014)	2.4E-10	0.28(0.058)	8.9E-06	0.69(0.080)	2.1E-14	−0.09(0.082)	0.28	9.4	4.8E-21	+ +
GNB1	1	0.05(0.009)	2.0E-08	0.23(0.059)	1.8E-04	0.42(0.085)	1.6E-06	0.06(0.084)	0.51	7.3	3.0E-13	+ +
KLF13	15	0.10(0.017)	1.1E-08	−0.01(0.060)	0.94	0.01(0.086)	0.89	−0.02(0.084)	0.80	4.9	8.5E-07	+ +
MYL5	4	0.09(0.017)	2.3E-08	0.20(0.059)	1.3E-03	0.45(0.083)	1.3E-07	−0.04(0.083)	0.60	7.5	5.0E-14	+ +
NINJ2	12	0.08(0.013)	4.3E-09	0.14(0.060)	0.03	0.24(0.087)	0.01	0.05(0.085)	0.59	6.3	2.2E-10	+ +
PRMT2	21	0.06(0.010)	3.6E-09	0.18(0.060)	0.01	0.27(0.087)	6.7E-03	0.09(0.085)	0.33	6.4	1.2E-10	+ +
SLC7A10	19	−0.27(0.042)	1.4E-10	−0.21(0.057)	7.4E-04	−0.31(0.082)	3.3E-04	−0.11(0.081)	0.18	−7.3	2.0E-13	− −
TPMT	6	0.10(0.013)	8.0E-15	−0.04(0.060)	0.49	−0.03(0.087)	0.78	−0.06(0.084)	0.49	6.4	1.1E-10	+ −

The trans-effects explain a substantial portion of the genetically-regulated variation in GWST expression levels. Our heritability estimates of GWST-gene expression levels ranged from 0.13 to 0.79: the rs4731702 trans-eQTL explained between 3-7.8% of the variance in expression, corresponding to 6-25% of the heritability (Table 2). Expression levels of the ten GWST genes are moderately-correlated in adipose tissue, with a mean ∣pairwise rho∣ of 0.29 (stdev = 0.15). SLC7A10 is the only GWST gene down-regulated by the T2D-risk allele (and hence the only transcript showing anti-correlated expression levels within the GWST genes). This pattern is consistent with the known ability of the KLF family of transcription factors to act as both transcription activators and repressors6.

Table 2

Heritability and trans-eQTL variance of GWST gene expression

Probe ID	Gene	CHR	Transcriptionstart site(b36)	%varianceexplainedin trans	h2	%h2 explainedin trans
ILMN_1767816	APH1B	15	61356801	3.1%	0.13	23.2%
ILMN_1684873	ARSD	X	2832010	6.6%	0.52	12.6%
ILMN_1693862	C8orf82	8	145722410	6.9%	0.27	25.4%
ILMN_1760320	GNB1	1	1706588	4.1%	0.19	21.3%
ILMN_1679929	KLF13	15	29406374	4.3%	0.70	6.2%
ILMN_1746948	MYL5	4	661710	4.0%	0.40	10.0%
ILMN_1731745	NINJ2	12	543722	4.7%	0.48	9.8%
ILMN_1675038	PRMT2	21	46879954	4.9%	0.55	8.9%
ILMN_1681087	SLC7A10	19	38391409	5.4%	0.79	6.8%
ILMN_1740185	TPMT	6	18236523	7.8%	0.48	16.4%

Further support for the hypothesis that the trans-effects are mediated by KLF14 expression comes from analysis of transcription-factor binding-sites in trans-associated genes using PSCAN7 with the JASPAR database8. KLF14 itself is not represented in JASPAR, but other KLF family members have closely related binding sites (and in some cases have been shown to compete for the same binding site)9, and KLF4 (the only KLF family member in the JASPAR database) and KLF14 share highly similar DNA binding C-terminal regions10. Though we found no evidence for enrichment after correction for multiple-testing when examining the 10 GWST genes alone, inclusion of a larger number of trans-associated genes (46 with trans p < 10−4 or 121 with trans p < 10−3) revealed strong evidence of enrichment for KLF4 binding sites. KLF4 was the most over-represented binding site in the former set (Bonferroni-corrected p = 0.01) and the second most over-represented site in the latter set (Bonferroni corrected p = 1.3 ×10−7) after EGR1. These data indicate that one feature of the transcripts showing trans-associations with the KLF14 SNPs is enrichment for KLF binding sites.

Having demonstrated that the same set of SNPs influences cis-expression of KLF14, trans-expression of members of the GWST-gene network, and a variety of metabolic traits including T2D and HDL-cholesterol, we sought to clarify the causal connections between these effects, and in particular to establish whether or not the trans-effects were likely to be mediating the metabolic associations at KLF14. First, we examined the correlations between trans-gene expression and concurrently-measured metabolic phenotypes. At an array-wide Bonferroni threshold of p< 1.9×10−6, expression levels of six of the ten GWST genes are associated with BMI and HDL-cholesterol, five each with triglycerides and fasting insulin levels, four with HOMA-IR (an index of insulin sensitivity) and two each with fasting glucose and adiponectin (Table 3). Compared to all genes on the array, this represents an enrichment for expression/metabolic phenotype associations, with significance ranging from p = 0.001 to p = 3.3 ×10−5. The strength of these associations is consistent with a causal link between trans gene expression and metabolic phenotypes, and provides clues to the biological processes in which these genes may participate.

Next we examined large-scale association data made available by trait-specific GWAS meta-analysis consortia, focusing on SNPs in the 250kb surrounding each GWST gene. The rs4731702 T2D risk-allele is associated with higher fasting insulin1, indicating that the primary effect on diabetes-risk is mediated by decreased peripheral insulin sensitivity. Accordingly, we focused on a set of insulin-resistance related traits including fasting insulin11, fasting glucose11, HOMA-IR11, T2D1, lipids (HDL, LDL, triglycerides)2, body fat distribution (BMI-adjusted WHR)12 and BMI13. In GWAS datasets ranging in size from 22,044-123,865 individuals, we found eight associations in five genes at a study-wide significance threshold of 1.03 ×10−4 (Table 4). (See methods for threshold determination). For example, SNPs near APH1B are associated with HDL (rs2729787; p=9.8 ×10−9) and triglycerides (rs17184382; p=1.5 ×10−5), and SNPs near KLF13 with BMI-adjusted WHR (rs4779526; p=1.8 ×10−5) and LDL (rs8034505; p=5.8 ×10−5). In addition, SNPs in MSRA (expression levels of which marginally failed to reach genome-wide significance: trans-association p=5.1 ×10−8) have been previously associated with waist circumference10, and are here associated with triglycerides (rs615171; p=7.5 ×10−7). This pattern of association signals reveals that variation involving GWST-genes has the potential to impact on insulin-resistance related traits, and thereby supports the notion that a subset of these genes are directly implicated in mediating the effects of KLF14 variation on disease-susceptibility.

Table 4

GWA meta-analysis signals (p < 1.03× 10−4) within 250KB of genome-wide significant trans genes. This table also includes the results for MSRA for which the trans association marginally failed to reach genome-wide significance (p = 5.1×10−8).

Gene	Trait	SNP	EffectAllele	Zscore	Pvalue
APH1B	HDL	rs2729787	T	5.73	9.81E-09
APH1B	Triglycerides	rs17184382	A	4.32	1.58E-05
C8orf82	Type 2 Diabetes	rs2294120	A	1.14 (Oddsratio)	8.43E-05
NINJ2	LDL	rs2302408	T	3.91	9.06E-05
SLC7A10	HDL	rs8182584	T	−5.11	3.19E-07
SLC7A10	Waist-Hip Ratio	rs7251505	A	−4.85	3.21E-06
KLF13	LDL	rs8034505	A	4.02	5.85E-05
KLF13	Waist-Hip Ratio	rs4779526	A	4.50	1.79E-05
MSRA	Triglycerides	rs615171	T	−4.95	7.50E-07
MSRA	Waist Circumference	rs7826222	G	5.75	8.89E-09

One of the more interesting transcripts revealed by these analyses is SLC7A10, a member of the solute carrier family that mediates transport of neutral amino acids. Adipose expression of SLC7A10 is highly heritable (h2=0.79) and is down-regulated by the KLF14 T2D risk-allele. SLC7A10 expression is strongly-associated with diverse metabolic phenotypes; negatively correlated with BMI (p = 3×10−48), insulin (p=1.1×10−51), HOMA-IR (p=7 ×10−48), glucose (p=6×10−7), and triglycerides (p =1 ×10−34) and positively correlated with HDL (p=7×10−30) and adiponectin (p= 1 ×10−12) levels (Table 3). The SLC7A10 locus contains independent (r2=0.03) SNPs associated with HDL (rs8182584; p=3.2×10−7), and BMI-adjusted WHR (rs7251505; p=3.2×10−6). The former SNP (rs8182584) is weakly associated to insulin (p=0.002) and BMI (p =1.4 ×10−3) suggesting that this gene has a wide-ranging role in metabolism.

Table 3

Association between expression of GWST genes and concurrently measured metabolic phenotypes. Values in each cell represent p value and (beta value).

Gene	Adiponectin	HDL	LDL	Trigylcerides	BMI	FastingInsulin	FastingGlucose	HOMA-IR
APH1B	7.7E-07 (−0.02)	1.9E-08 (−0.11)	0.01 (0.02)	5.2E-17 (0.13)	5.2E-15 (0.01)	6.7E-18 (0.002)	0.001 (0.05)	1.2E-12 (0.05)
ARSD	0.02 (−0.008)	6.8E-09 (−0.10)	0.02 (0.02)	3.6E-11 (0.09)	5.3E-11 (0.01)	4.3E-16 (0.002)	4.2E-10 (0.07)	2.1E-11 (0.04)
C8orf82	2.8E-05 (0.02)	4.4E-08 (0.11)	0.008 (−0.03)	4.1E-05 (−0.07)	2.4E-14 (−0.01)	6.0E-05 (−0.001)	0.52 (−0.009)	0.0003 (−0.02)
GNB1	0.004 (−0.008)	6.1E-10 (−0.08)	0.03 (0.01)	2.1E-12 (0.08)	2.9E-21 (0.01)	8.4E-08 (0.001)	0.00053 (0.03)	9.8E-06 (0.02)
KLF13	0.17 (0.005)	0.11 (−0.04)	0.86 (−0.002)	0.33 (0.02)	0.93 (0)	0.87 (0)	0.02 (0.04)	0.48 (−0.0002)
MYL5	0.00075 (0.02)	0.00053 (0.10)	0.19 (−0.02)	0.8 (−0.005)	8.4E-05 (−0.009)	0.007 (−0.001)	0.47 (−0.01)	0.001 (−0.02)
NINJ2	0.85 (−0.002)	0.03 (−0.04)	0.03 (0.02)	0.00012 (0.06)	0.003 (0.005)	0.0006 (0.001)	0.51 (0.009)	0.05 (0.02)
PRMT2	0.04 (0.005)	0.06 (0.03)	0.56 (−0.004)	0.21 (−0.015)	1.1E-05 (−0.006)	0.008 (−0.0004)	0.3 (−0.01)	0.002 (−0.01)
SLC7A10	1.6E-14 (0.07)	7.8E-30 (0.65)	2.7E-05 (−0.12)	1.7E-34 (−0.55)	3.6E-48 (−0.08)	1.2E-51 (−0.01)	6.7E-07 (−0.21)	5.3E-45 (−0.27)
TPMT	0.0003 (−0.01)	2.3E-09 (−0.11)	0.01 (0.023)	4.9E-16 (0.12)	7.5E-24 (0.02)	7.3E-21 (0.002)	0.00057 (0.05)	6.7E-15 (0.06)

Our data provide convincing evidence of a bona fide adipose trans-eQTL and implicate this trans-expression network in the link between KLF14 variation and risk of metabolic disease. The trans-regulation uncovers novel biological links between previously-identified genome-wide significant associations at KLF14 (HDL; T2D), APH1B (HDL) and MRSA (Waist circumference) and to additional signals where metabolic trait associations have not yet been established to genome-wide significance (SLC7A10, KLF13, C8orf82, NINJ2). These links provide a framework for hypothesis-directed investigation of genetic interactions among GWAS loci and provide an example of the power of ‘integrative genomics’ to leverage ‘omics data from multiple sources to discover new biological and functional insights.