Expression quantitative trait analyses to identify causal genetic variants for type 2 diabetes susceptibility

Abstract

Type 2 diabetes (T2D) is a common metabolic disorder which is caused by multiple genetic perturbations affecting different biological pathways. Identifying genetic factors modulating the susceptibility of this complex heterogeneous metabolic phenotype in different ethnic and racial groups remains challenging. Despite recent success, the functional role of the T2D susceptibility variants implicated by genome-wide association studies (GWAS) remains largely unknown. Genetic dissection of transcript abundance or expression quantitative trait (eQTL) analysis unravels the genomic architecture of regulatory variants. Availability of eQTL information from tissues relevant for glucose homeostasis in humans opens a new avenue to prioritize GWAS-implicated variants that may be involved in triggering a causal chain of events leading to T2D. In this article, we review the progress made in the field of eQTL research and knowledge gained from those studies in understanding transcription regulatory mechanisms in human subjects. We highlight several novel approaches that can integrate eQTL analysis with multiple layers of biological information to identify ethnic-specific causal variants and gene-environment interactions relevant to T2D pathogenesis. Finally, we discuss how the eQTL analysis mediated search for “missing heritability” may lead us to novel biological and molecular mechanisms involved in susceptibility to T2D.

Key Words: Type 2 diabetes; Single nucleotide polymorphisms; Expression quantitative trait locus; Expression regulatory SNPs; Gene-environment interaction; Genome-wide association study

Core tip: Identification of genetic variants that modulate the susceptibility to disease and elucidating their function at the molecular level is a major focus of type 2 diabetes (T2D) research. This article highlights the utility of expression quantitative trait analysis in discovering regulatory variants that increase susceptibility to T2D by modulating the expression of transcripts in tissues important for glucose homeostasis.

GENETIC DISSECTION OF TYPE 2 DIABETES SUSCEPTIBILITY

Diabetes is one of the most prevalent metabolic disorders, characterized by elevated levels of plasma glucose, and is responsible for significant mortality and morbidity in human populations worldwide[1]. The latest estimate from the International Diabetes Federation indicates a global prevalence rate of 8.4% in adults and 382 million cases of diabetes in 2013[2]. It is one of the common diseases with a well-accepted genetic contribution[3]. Type 2 diabetes (T2D), a late onset subtype of diabetes, results from a derangement in the complex interplay of multiple physiological processes known to be involved in systemic glucose homeostasis. These processes include peripheral glucose uptake in muscle, secretion of hormones and incretins form pancreas and intestine, secretion of cytokines/adipokines from adipose tissue, hepatic glucose production, and neuro-endocrine regulation by central nervous system[4,5]. However, the relative contribution of these processes to T2D pathogenesis is debated. Based on this knowledge on intertwined and complex physiological processes it can be anticipated that T2D is a heterogeneous conglomeration of phenotypes, caused by multiple genetic perturbations and affecting different biological pathways. Predictably, deciphering the genetic etiology of T2D has remained challenging.

Until the last decade, searching for an association between T2D and sequence variants of selected candidate genes was the mainstay of research for finding genetic susceptibility factors. Based on available technology in those studies, researchers selected candidate genes either from loci detected by genome-wide linkage analyses or based on known physiological functions. In our earlier reviews, we discussed the knowledge gained from such studies in detail[6,7]. Success from those endeavors was very limited. However, this approach has identified genetic variants in the TCF7L2 gene, to date is the best replicated and strongest (relative risk approximately 1.4) genetic susceptibility factor for T2D[8], but its role is still controversial[9-11].

In the middle of the last decade, a transformative change took place in the field of genetics of complex disease research. Advances in high-throughput genotyping technology, availability of the complete human genome sequence, a dense catalogue of common genetic variants, and a population-specific linkage disequilibrium map of these variants lead to the implementation of genome-wide association studies (GWAS), which interrogate the entire genome to identify common genetic variants (minor allele frequency ≥ 0.05) associated with a disease[12]. GWAS have yielded unprecedented success in identifying well-replicated susceptibility loci for T2D, glucose homeostasis traits, obesity, and related metabolic phenotypes[3,13-15]. Nevertheless, these successes come with significant caveats. Based on the most recent analyses, the 63 T2DM-associated loci discovered so far in Caucasian populations together account only for 5.7% of the liability-scale variance in disease susceptibility, and sibling relative risk (λ_s) attributed jointly by these variants is 1.104[13]. Moreover, few of the T2D loci identified primarily in European- or Asian-derived populations are convincingly replicated in African American, Native American, and Hispanic populations, all of whom have a higher prevalence of T2D than Caucasians[14,16]. These GWAS-identified loci do not appear to explain the well-established roles for adipose, muscle, and liver in diabetes pathogenesis[17], and few of these loci have been linked to a molecular mechanism. Several investigators have attempted to implicate function to T2D-associated loci based on their proximity to a gene, assuming that the associated single nucleotide polymorphisms (SNP) alters the function of a nearby gene[18]. Some have drawn enthusiastic conclusions about the role of these variants exclusively in insulin secretion[19]. However, proof of such an assumption is lacking. Given the small effect on T2D susceptibility and the statistical noise inherent in performing 10⁶ or more tests, exclusive reliance on larger T2D GWAS alone is unlikely to identify the source of undefined T2D susceptibility (often referred to as “missing heritability”[20]).

EXPRESSION QUANTITATIVE TRAITS: MOLECULAR ENDOPHENOTYPES

One of the major findings from the T2D GWAS is that most of the trait-associated SNPs are located in intronic, intergenic, or other non-coding regions of the genome[3,21]. Further fine mapping analysis also failed to find any coding or other variants that would provide a molecular biological explanation of the elevated disease risk attributed by these loci.

The abundance of a transcript is a quantitative trait. Studies in human populations showed a wide, heritable variation of transcript levels among individuals, and thus lead to the concept of “expression quantitative trait loci”(eQTL)[22,23]. The heritability of eQTLs has been replicated in multiple human tissue or cell types, with approximately 30% of eQTLs having h² > 0.3, and an estimated 58%-85% being heritable[24-28]. The abundance of a transcript can be directly modified by polymorphisms in non-coding regulatory elements. Many SNPs are associated with quantitative transcript levels and are considered as expression regulatory SNPs (eSNPs). eSNPs close to the transcription start sites (TSS) of the eQTLs are named “cis” or “local” eSNPs , whereas eSNPs located > ± 500 kb from the TSS or on a different chromosome are considered “trans” or “distal” eSNPs[22,29]. Similar to a published study[30], here we will refer to eQTLs as the transcripts rather than SNP-transcript pair, and eSNPs as the genetic variants (SNPs) associated with the expression profile of a transcript.

Based on this knowledge, many laboratories (including ours) hypothesized that GWAS-associated non-coding variants are eSNPs and can modulate T2D susceptibility by altering transcript levels (or splicing). This concept is based on the “central dogma” of gene expression and presents a causal model of genetic susceptibility (Figure 1). In this model, transcript abundance is considered as an intermediate phenotype between genetic loci (DNA sequence variants) and subclinical (e.g., insulin resistance) or clinical (e.g., T2D) phenotypes. Since transcript abundance is a proximal molecular endophenotype affected by genetic variants, it is likely to be a less heterogeneous phenotype (compared to complex clinical phenotypes like those of T2D), and thus more amenable to genetic mapping methods due to superior statistical power.

Figure 1 A causal model of genetic susceptibility. Genetic regulatory architecture modulates molecular phenotypes in interaction with environmental factors and alters disease susceptibility. eSNP: Expression regulatory SNP; eQTL: Expression quantitative trait loci; T2D: Type 2 diabetes.

EQTL MAPPING

Study designs and analytical frameworks for eQTL mapping are similar to those for mapping any other quantitative traits [e.g., body mass index (BMI), fastin glucose, glycosylated hemoglobin]. However, genetic analysis of human phenotypes including QTLs carries a unique set of problems[29]. In general, eQTL analyses integrate genome-wide expression (in tissues or cells) and genotype data in multiple individuals (related or unrelated). These analyses use linkage- or association-based statistical genetic methods to map regulatory regions and genetic variants that may explain individual variations in transcript expression. Microarray- or RNA-seq[31-33] based methods are used to generate large numbers of quantitative transcript phenotypes. Therefore, the number of statistical tests involved in eQTL mapping studies is significantly higher than in traditional QTL analysis[34]. A detailed discussion on methods used in eQTL analysis is beyond the scope of this article, and we refer our readers to other reviews on this specific subject[29,34-36].

Published eQTL studies have implemented linkage analysis by using 400-2000 microsatellite makers[24,26] to localize regulatory intervals, whereas other studies have genotyped large numbers of common SNPs ( > 100000) to discover the eSNPs[25,28,37] associated with eQTLs. With the advancement of genomic technology, we can now simultaneously genotype more than 4.5 million SNPs or can have a whole genome sequence for each individual included in an eQTL study by highly multiplexed “next generation” sequencing methods[38]. These advances pose additional statistical and computational challenges, and will require appropriate correction and adjustment of significance thresholds for the massive number of independent tests performed (and hypotheses tested) to control false discovery. The power to detect eSNPs depends on their effects (average difference in the transcript abundance between genotypes, scaled by the standard deviation of the transcript abundance within genotype classes) and allele frequency[34]. Consequently, detection of eSNPs with a lower effect allele frequency and a lower effect size will require a larger sample size.

One interesting observation from published eQTL studies is that most of the strong eSNPs are located near the TSS with no discernable trend in the 5’ or 3’ direction[28,39,40]. As a result, most studies consider SNPs within close proximity of the TSS (± 500 kb window) as cis-eSNPs. Since the genomic context of most eQTL transcripts are known, statistical adjustment for the actual number of SNPs tested within 500 kb will be more appropriate for cis-eSNP discovery. Any SNP outside the cis-region is tested as a trans-eSNP for a transcript. The molecular biological basis of trans-regulation is less studied; current information suggests that the variants that affect transcription factors, miRNAs, or long-range chromatin interaction may act as trans-eSNPs. To identify trans-eSNPs, the number of tests needed is far greater, and the tests require more stringent significance threshold criteria and a larger sample size. Thus, use of a false discovery rate based on a permutation analysis to correct for multiple testing[34], and considering the correlation among transcript levels and highly correlated SNP structures, are useful approaches to identify this biologically important class of regulatory SNPs.

Several heterogeneous sources of variability hidden in the data may lead to both spurious eSNPs and missed associations in eQTL analyses if not properly addressed. Statistical models that correct for hidden structures within the sample (such as race, admixture, and family relatedness), artifacts in expression data (including batch effects and probe bias), environmental influences, and other known and unknown factors are required to improve sensitivity and interpretability of eQTL analyses[41]. Methods that showed significant usefulness in tackling these confounding factors include Bayesian approaches developed by Stegle et al[42] (implemented in probabilistic estimation of expression residuals or probabilistic estimation of expression residuals software), linear mixed-effects model-based approaches developed by Listgarten et al[43] (implemented in LMM-EH-PS or Linear Mixed Model-Expression Heterogeneity-Population Structure software), surrogate variable analysis, and inter-sample correlation emended approaches[44,45].

The heavy computational burden involved in eQTL analyses sometimes forces researchers to restrict their analysis to a small subset of selected transcripts and SNPs. Improvement of computational algorithms, parallelization of programs by efficient scripting, and utilization of efficient processing hardware are among many approaches needed to improve scalability and computational efficiency required for eQTL analyses. Implementation of these approaches will enhance discovery by increasing the capacity to utilize the complete data set[46,47].

EQTLS AND DISEASE GENE MAPPING

Molecular and cell biological experiments in model organisms and cells have significantly advanced our understanding about the role of non-coding DNA sequences in genetic regulation, transcriptional circuitry, the transcriptional apparatus, and chromatin regulation. This work has led to new insights into the complex mechanisms involved in dysregulation of gene expression in various human diseases[48]. Recent genome-wide studies in human cells by different international consortia [including ENCyclopedia Of DNA elements (ENCODE)][49] further have improved our mechanistic understanding of the role of DNA sequence variants in quantitative modulation of gene expression[50-52]. eQTL studies have been extensively used to identify genetic regulators involved in natural variation of gene expression[28,37,39] and to understand tissue-specific architecture of genetic regulatory mechanisms[24,30,53-59].

However, an intriguing application of eQTL mapping is the use of eSNP data to interpret disease or disease-related phenotypic association signals, and thereby elucidate specific biological mechanisms underlying the increased genetic risk attributed by the DNA sequence variants. Identification of genetic variants simultaneously associated with disease and eQTLs (in relevant tissue) significantly facilitates identification of potential causal genes. Discovery of genetic variants in ORMDL3 as a susceptibility factor for childhood-onset asthma[60] and VNN1 variants that influence high-density lipoprotein cholesterol concentrations[26] are two early examples of the successful implementation of eQTL mapping in disease gene hunting. The review by Cookson et al[36] offer a more detailed discussion on those success stories.

Several recent studies have integrated GWAS and eQTL analyses (data generated in different sets of subjects) and have used the overlap of two signals as a tool to interpret GWAS findings. Although this work is a good starting point, we need to be cautious about using the overlap of two statistical signals (eSNP and the disease phenotype-associated SNP/phSNP). Careful thought is required before making a claim of identifying a disease-causing variant. Montgomery and Dermitzakis (2011) described three situations[41] when a coincidence of eQTL and disease phenotype GWA signal may distract from identification of causal variants: (1) eSNP and phSNP are in the same linkage disequilibrium (LD) block but are two different SNPs. This is not considered as exact overlap, and they may tag different causal variants; (2) eSNPs and phSNPs are the same but SNP density differs between the eQTL and GWAS data. Lack of proper resolution in one or both studies may be misleading and will not elucidate the correct functional SNP; and (3) eSNPs may have a pleiotropic effect and may regulate the expression of “gene Y” in “tissue 1”, but the same eSNP may regulate the expression of “gene X” in “tissue 2”. Thus, if the eQTL study is done in “tissue 1” (a “surrogate” tissue) but not in “tissue 2” (the “disease-relevant” tissue in which the true causal effect is manifested), then despite the overlap of eSNPs and phSNPs, we will incorrectly link “gene Y” to the disease phenotype.

In general, eSNPs that are universal have a stronger effect, but a significant proportion of eSNPs show tissue-specific effects[30,53,54]. However, it is difficult to select “relevant” tissue, or the relevant tissue may not be accessible from human subjects for analysis for many complex diseases. Ongoing efforts of international consortia, including GTEx, to develop multi-tissue eQTL databases (Table 1) is a significant step forward in addressing this limitation[61-64].

Table 1 Selected expression quantitative trait loci databases.

Database	Website (URL)	Cell/tissue type	Project
eQTL Browser	http://eqtl.uchicago.edu/cgi-bin/gbrowse/eqtl/	LCL, liver, brain, fibroblast, T-cell	17 projects
Genvar	http://www.sanger.ac.uk/resources/software/genevar/	Adipose, LCL Skin fibroblast from healthy female twins	MuTHER
		LCL from 8 populations	Hapamap3
		Fibroblast, LCL and T-cell from umbilical cord	GenCord
GTEx eQTL Browser	http://www.ncbi.nlm.nih.gov/gtex/GTEX2/gtex.cgi	Multiple tissues including liver, brain regions, LCL	GTEx
PACdb	http://www.pacdb.org/	Gene-drug or GXD eSNPs from LCL model	Dolan and Cox lab Lusis lab
SGR Database	http://systems.genetics.ucla.edu/	22 mouse and several human datasets.
		Data includes aortic endothelial and smooth muscle, adipose, brain, liver, macrophages and muscle tissue
		Includes GXE eSNP data from cell experiments
SCAN	http://www.scandb.org/newinterface/about.html	CEU and YRI LCLs from HapMap	Cox Lab
seeQTL	http://www.bios.unc.edu/research/genomic_software/seeQTL/	HapMap LCLs

SGR: Systems genetics resource; eQTL: Expression quantitative trait loci; T2D: Type 2 diabetes; eSNP: Expression regulatory SNP; LCL: Lymphoblastoid cell lines; GXD: Gene-by-drug interaction; GXE: Gene-by-environment interaction; CEU: HapMap caucasian from CEPH collection; YRI: HapMap African from Yuroba, Nigeria.

Many investigators have developed statistical approaches to formally test the overlap of GWAS and eQTL signals to distinguish accidental colocalization from true sharing of causal variants. The regulatory trait concordance method designed by Nica et al[65] accounts for local LD structure and integrates eQTL and GWAS results to reveal the subset of association signals due to cis- or trans-eQTLs. He et al[66] (2013) developed an algorithm named “Sherlock” based on a Bayesian statistical framework to identify potential gene-disease associations by matching genetic signatures of expression (collective information of cis- and trans-eSNPs) of a gene to that of the disease phenotype by using GWAS data of the disease and the eQTL data of related tissue. These novel approaches are likely to expand our ability to harvest new insights from genetic association studies for disease phenotypes.

T2D-ASSOCIATED VARIANTS ARE ESNPS IN TISSUES IMPORTANT FOR GLUCOSE HOMEOSTASIS

Genome-wide eQTL analyses in transformed lymphocytes (lymphoblastoid cell lines, or LCLs) provided the first evidence that SNPs associated with complex diseases phenotypes are more likely to be eSNPs than minor allele frequency-matched SNPs randomly selected from high-throughput GWAS genotyping platforms. Nicolae et al[67] (2010) utilized an Affymetrix GeneChip Human exome 1.0 ST array to generate exon-level expression data of LCLs from 87 Caucasian (CEU) and 89 African (YRI) subjects from the HapMap project. They performed a quantitative-trait transmission disequilibrium test to identify eSNPs from 2 million genotyped SNPs. A study by Nica et al[65] (2010) utilized an Illumina Sentrix WG-6-V2 whole-genome expression array to generate total transcript-level expression data of LCLs from 109 unrelated CEU subjects (from the HapMap 3 project) and performed Spearman rank correlation analysis to identify eSNPs from 1186075 genotyped SNPs. Key findings from these studies[65] include: (1) SNPs reproducibly associated with complex human traits are likely to be eSNPs; (2) Enrichment of complex trait GWAS-implicated SNPs are more evident among cis-eSNPs but not among trans-eSNPs; and (3) eSNPs discovered in LCLs are more strongly enriched for SNPs associated with immunity-related conditions (e.g., Crohn’s disease, type 1 diabetes, rheumatoid arthritis), but such enrichment was not observed for metabolic disorders (e.g., T2D and coronary artery disease). These studies indicate that eQTL studies using surrogate tissue samples may be helpful for some diseases. However, understanding the functional role of T2D-associated SNPs will probably require expansion of eQTL studies into tissues more relevant for T2D pathophysiology. These studies also had significantly lower power to identify trans-eQTLs due to comparatively small sample sizes, and will require reevaluation of the role of trans-eSNPs in larger sample sets.

Zhong et al[68] (2010) used genetics of gene expression (GGE) analysis in tissues from two cohorts of human subjects (Cohort 1: liver-specific GGE cohort with post mortem liver samples from 427 subjects; Cohort 2: liver, subcutaneous adipose and omental adipose from 922 subjects who had Roux-en-Y gastric bypass surgery). They identified 18785 unique eSNPs in the combined set of data. They found 2189, 2286, and 1999 eSNPs specific to liver, omental adipose, and subcutaneous adipose, respectively. However, they also noticed that 72% of cis-eSNPs identified in liver, 79% of those found in omental adipose and 80.5% from subcutaneous tissue were also found in the other two tissues. Given the metabolic relevance of these tissues, they further interrogated data from three large-scale T2D GWAS datasets to test whether the set of eSNPs were more likely to be associated with T2D compared to randomly selected SNPs. These tissue eSNPs showed a significant enrichment of T2D-associated SNPs. For example, in the DIAGRAM (DIAbetes Genetics Replication and Meta-analysis) GWAS data set, 7.34% of the eSNPs showed a significant association with T2D (P < 0.05) compared to an average of 6.12% SNPs in the random sets, representing a modest 1.20-fold enrichment for SNPs in the eSNP (or SNP in LD at r² > 0.89) set over the random sets (p-enrichment = 1.33 × 10^-9)[68]. In that study, omental adipose tissue eSNPs also showed further significant enrichment when restricted to adipose expression network genes differentially expressed with T2D. Thus, these studies support the notion that T2D- associated SNPs may modulate expression of transcripts in tissues relevant for glucose homeostasis.

Fu et al[53] (2012) analyzed eQTLs in blood (n = 1240) and other tissues (liver, n = 62; muscle, n = 62; subcutaneous adipose, n = 83; and visceral adipose, n = 77); out of 1954 SNPs associated with complex disease traits from a GWAS catalogue, 907 were cis-eSNPs. However, 28.7% of these trait-associated cis-eSNPs showed a tissue-specific (in blood versus other tissue) and discordant effect on gene expression. The discordant effect includes tissue-specific regulation, alternative regulation by different eSNPs, different effect size and, in a few cases, opposite allelic direction. The study also showed that SNPs associated with complex traits are more likely (P = 2.6 × 10^-10) to exert a tissue- specific effect on gene expression[53]. No comparisons were made between other tissues due to small sample size. This study indicates that use of tissues in eQTL analysis may have implications for inferring transcriptional effects of SNPs, especially for the complex disease susceptibility variants.

This work also emphasizes the importance of investigating disease-relevant tissue for characterizing functional effects of T2D and other disease-associated variants. However, it is difficult to determine “relevant tissue” even for diseases with known pathophysiology. T2D is clearly of polygenic etiology, and relevant tissue could be distinct for genes involved. Moreover, gene expression is regulated by environmental (e.g., diet), epigenetic, and other unknown factors, and eQTL discovery from tissue samples may be affected by the physiological state of the donors[41]. For example, profound hyperglycemia and dyslipidemia observed in T2D subjects will modulate and even may mask primary causal changes in genetic regulatory networks. Thus, multi-tissue eQTL analysis in physiologically characterized individuals could be a safe option to scrutinize the circularity of cause and effect in genetic regulatory signals, and holds the promise to offer insights into the novel mechanisms driving genetic susceptibility to T2D.

Most initial eQTL studies seeking to identify a regulatory role for T2D-associated SNPs have focused on cis-eQTLs. However, studies by Voight et al[69] (in adipose, n = 603; and blood, n = 745 subjects) and our laboratory (in adipose and muscle of 168 non-diabetic subjects who were physiologically evaluated) showed that only a few top T2D GWAS-identified signals can be explained as cis-eQTLs, and T2D-associated non-coding SNPs are less likely to regulate expression of the closest gene[70]. Results were similar in an eQTL analysis that used human islet cells from 63 cadaver donors[71]. A genome- wide study by our laboratory[72] in adipose and muscle tissue of 62 subjects (31 insulin- resistant and 31 insulin-sensitive subjects matched for BMI) showed that at a less stringent threshold (P < 0.0001), among 68 well-replicated T2D/glucose homeostasis-associated SNPs, 25 and 19 of them were eSNPs in adipose and muscle, respectively (Figure 2). However, after stringent (Bonferroni) correction, only SNP rs13081389 was a cis-eSNP for the SYN2 gene in adipose (P < 4.7 × 10^-8, 15507 expressed transcripts were tested in adipose). Interestingly, these 68 SNPs showed significant enrichment for trans-eSNPs in adipose and muscle, but not in LCLs[72]. Many of these trans-eSNPs show associations with expression of ≥ 10 transcripts and may be a “master regulator”. Expanding this search for the top 1000 T2D-associated SNPs from a Wellcome Trust Case Control study also confirmed the trans/distal regulatory SNPs[72]. We also showed that replicated T2D- and glucose homeostasis-associated SNPs are enriched for trans-eQTLs for transcripts differentially expressed between insulin-resistant and insulin-sensitive people[72]. A recent eQTL study using a large cohort of blood samples also supported the trans-regulatory role of 233 complex trait-associated SNPs[73]. Thus, the genetic regulatory architecture of T2D is complex, tissue-specific, and likely extends beyond the cis-regulatory mechanism.

Figure 2 Type 2 diabetes or glucose homeostasis traits associated variants are expression regulatory SNP. We tested Cis and Trans regulatory role of 68 SNPs that showed reproducible associations with T2D or Glucose homeostasis traits[72]. At a threshold of P < 0.0001, 25 and 19 of these SNPs in adipose and muscle, respectively, showed association with expression of a cis- or trans-transcript. This figure represents a CIRCOS plot of eQTL and eSNP chromosomal location relationships, indicating the predominance of trans-regulation among 183 and 62 significant (P < 0.0001) eQTL-eSNPs associations in adipose and muscle respectively. Rare cis-regulation (SYN2 in adipose and JAZF1 in muscle) is marked. eSNP: Expression regulatory SNP; eQTL: Expression quantitative trait loci; T2D: Type 2 diabetes.

EQTL ANALYSIS FOR PRIORITIZING T2D-ASSOCIATED VARIANTS TO IDENTIFY NOVEL CANDIDATE GENES

The multiple testing corrections utilized in genome-wide statistical analyses allow detection of only the strongest effects and penalize weaker associations that may be biologically meaningful[74]. Investigators have implemented several approaches to prioritize T2D association signals from large GWAS datasets to identify biological mechanisms responsible for genetic predisposition. One common approach includes selection of genes close to T2D GWAS-implicated SNPs and shows differential expression in T2D subjects compared to normoglycemic subjects (or in animal models of T2D). This approach is based on the idea that T2D-associated variants may modulate the expression of nearby genes in tissues important for glucose homeostasis. Parikh et al[75] used publicly available expression microarray data from different tissues (pancreas, adipose, muscle, and liver from T2D patients and rat models of T2D) to prioritize among the 275 genes located near 1170 T2D GWAS-implicated SNPs. A recent study by Taneera et al[71] used expression profiling of human pancreatic islet cells for functional prioritization of genes in the vicinity of 47 T2D-associated SNPs. However, available data from several human tissue eQTL analyses indicate that only a few T2D-associated SNP act as cis-eSNPs, and no enrichment of differentially expressed genes was observed around T2D GWAS-implicated variants[72]. Thus, a logical alternative for prioritizing T2D-associated variants is to utilize a reverse genetics approach and restrict the genetic search space to the subset of variants that are eSNPs in relevant tissues. These eSNPs are statistically associated with expression of transcript and thus have a strong possibility of being a “key driver” in perturbing gene-expression regulatory networks.

Selecting the genes based on eSNPs among those also associated with T2D in large GWAS datasets will prioritize genes with a significantly high chance of being causally involved with susceptibility to T2D, and thus may be helpful in identifying additional genetic susceptibility loci from GWAS datasets. A genome-wide analysis of adipose tissue transcriptomes from 62 insulin-resistant and -sensitive subjects identified 172 differentially expressed transcripts[76]. We checked adipose eQTL data from the MuTHER study[55] to find eSNPs of these differentially expressed transcripts. We further mined the DIAGRAM GWA meta-analysis results[13] for association of these eSNPs with T2D. This analysis[77] identified that the strongest cis-eSNP (rs11037579, P = 4.21 × 10^-6) for the HSD17B12 in adipose tissue was also associated with T2D [P = 3.80 × 10^-4, OR = 1.06 (95%CI: 1.03-1.1)]. Individuals carrying the T2D risk allele T for the intronic SNP rs11037579 had lower expression of HSD17B12 in adipose tissue. This result corroborates the finding that HSD17B12 expression is downregulated in the adipose tissue of insulin-resistant subjects. The HSD17B12 gene codes a bifunctional enzyme involved in the biosynthesis of estradiol and the elongation of very long chain fatty acids. Several variants within ± 500 kb of this gene are eSNPs (including a 3’UTR SNP rs1061810) in adipose, LCL, and other tissues, and show an association with T2D (although below the genome-wide threshold) (Figure 3). Further functional studies will be required to identify true causal SNPs. However, this integrative approach demonstrates the validity of such an approach in prioritizing novel T2D susceptibility loci. In fact, two recent integrative genomic studies showed that eSNPs for PFKM (SNP rs11168327) gene in muscle and ARAP1 (SNP rs11603334) gene in pancreatic beta cell are associated with T2D[78,79].

Figure 3 Prioritizing type 2 diabetes-associated variants by expression quantitative trait loci analysis: An example. HSD17B12 is one of 172 genes differentially expressed in adipose tissue of insulin-resistant (IR, n = 31) vs insulin-sensitive (IS, n = 31) subjects in a genome-wide study (A) by Elbein et al[76]. Its expression in subcutaneous adipose of non-diabetic subjects (n = 141) also shows a significant correlation (B) with insulin sensitivity (SI). Strongest cis-eSNP for adipose tissue (C) expression of HSD17B12 (in adipose eQTL from the MuTHER project)[55] is also associated with T2D (D) in a large GWAS meta-analysis (in DIAGRAM.v3 data from 12171 T2D and 56862 controls)[13]. This locus also includes a 3’UTR SNP rs1061810 that shows association (E) with T2D and expression of HSD17B12 (in qRT-PCR analysis in adipose tissue from 141 non-diabetic subjects). eSNP: Expression regulatory SNP; eQTL: Expression quantitative trait loci; T2D: Type 2 diabetes.

EQTL AND BIOLOGICAL NETWORK ANALYSIS TO IDENTIFY ETHNIC-SPECIFIC GENES FOR T2D:

Age-standardized prevalence of T2D varies among ethnic and racial groups[14,80]. T2D is almost twice as prevalent in adult non-Hispanic African Americans (14.9%) in the United States compared to European Americans (7.6%)[81]. Yet only a few of the associated T2D-loci - identified primarily in European- or Asian-derived populations - are replicated in African American, Hispanics, and Native Americans[14,16,82-84]. Intriguingly, studies have identified distinctive physiologic features of glucose homeostasis in African Americans and Hispanics[85-87]. Compared to non-Hispanic Caucasians matched on age, gender, and BMI, African Americans are more insulin-resistant (lower S_I), but show a greater acute insulin response to intravenous glucose (AIR_G) and a higher disposition index (DI = S_I× AIR_G). A genetic basis for these physiological differences seems likely, but remains unidentified.

Published studies of expression across ethnic groups (mostly restricted to lymphocytes or HapMap LCLs) showed distinct ethnic-specific expression[37,88-90]. Zhang et al[90] (2008) reported differential expression of up to 67% of transcripts between LCLs from subjects of European (CEU) and African (YRI) descent, with enrichment of ribosome biogenesis, antimicrobial response and cell-cell adhesion. Spielman et al[37] (2007) attributed the 1097 genes that differed between CEU and Asian (CHB) LCL samples to eSNP frequency. Our comparison of genome-wide expression profiles (using an Agilent 44K expression array) from adipose and muscle tissue of non-diabetic Caucasians (n = 40) and African Americans (n = 22) identified transcripts associated with insulin sensitivity (S_I), many of which (e.g., CLIC6, HSD11B1, SERPINA3, THBS1, TMEM135 and TNMD in Adipose ) show distinct ethnic-specific expression[76].

Comparison of adipose tissue expression data between Caucasians and African Americans in a larger cohort (using an Illumina –HT12.V4 array for 99 Caucasians and 37 African Americans) identified 117 differentially expressed (fold change ≥ 1.5 and false discovery rate ≤ 5%) transcripts[91]. By mining adipose tissue eQTL data from the MuTHER project[55], we found that about 35% of these differentially expressed transcripts are strongly modulated (P < 1 × 10^-5) by cis-eSNPs in adipose tissue. In line with the findings by Spielman et al[37] (2007) in LCL, we also found that in adipose tissue, the degree of differential expression (fold change African Americans/Caucasians) shows strong concordance with the difference in the effect allele frequency of top cis-eSNPs (Figure 4) between HapMap African (YRI) and Caucasian (CEU) subjects.

Figure 4 Population differences in expression of transcripts in adipose tissue is accounted for by the effect allele frequency difference of expression regulatory SNPs among racial groups. X axis: Fold change in average expression of 41 transcripts between African-American (AA, n = 37) and Caucasian (CA, n = 99) subjects. Y axis: Differences in strongest eSNP allele frequency of these transcripts between HapMap subjects of Caucasian (CEU) and African (YRI) ancestry for alleles associated with higher expression. eSNP: Expression regulatory SNP.

These studies suggest that the distinct genetic architecture of eSNPs determines the ethnic-specific expression profile in tissues important for glucose homeostasis. Ethnic-specific derangements of gene expression networks in tissues involved in glucose homeostasis may explain distinctive physiologic effects, including differences in insulin action and secretion between ethnic and racial groups. Perturbation of gene expression networks associated with early pathophysiologic events (including insulin resistance) is driven by regulatory variants (eSNPs). The distinct genetic architecture of these variants (including linkage disequilibrium and allele frequency) may determine their ethnic-specific (or predominant) effect on expression and T2D susceptibility. Thus, integration of genome-wide expression analysis and eQTL analyses may be a useful approach to identify the primary genetic factors for ethnic-specific susceptibility to T2D.

Expression of transcripts involved in the same biological function tend to be co-regulated by similar factors (genetic or environmental) and can be identified as distinct network modules, where genes within a module are more highly interconnected (correlated) with each other than genes in other modules. Statistical approaches like weighted gene co-expression network analysis (WGCNA software package developed in “R” programming environment implements this analytical method) are useful for identifying modular structures of the co-expression networks[92,93] in tissues important for glucose homeostasis. Evaluation of the correlation of each module eigengene with the S_I and other T2D-related metabolic phenotypes, and determination of the preservation of these modules between ethnic groups based on observed network density and connectivity, will identify molecular processes or molecular interaction structures associated with phenotypes that undergo ethnic-specific reconfiguration by genetic or non-genetic causal regulators.

Several recently developed statistical metrics[94,95], including modular differential connectivity, offer powerful tools to identify the modules with significant ethnic-specific changes in interaction strength. The eSNPs are causal variants (or in linkage disequilibrium with causal variants) that regulate the expression level of neighboring (or distal) genes. Thus, eSNPs serve as a primary source of natural perturbation to infer causal relationships among and between genes in gene-expression networks[96]. The distinct allelic architecture of these SNPs may determine ethnic-specific modular differential connectivity. Genes with eSNPs can be considered as “parent nodes” in expression networks. This information is used as a “structure prior” in the network reconstruction analysis to orient the edges of the networks. Reconstructing ethnic-specific networks by utilizing different causality modeling methods, including Bayesian network reconstruction approaches, may identify key causal regulators of these networks[97,98]. Thus, a multiscale biological network analysis that utilizes eQTL information to distinguish causal from correlated disease effects is a novel approach to understand how causal regulators propagate their effects in mediating ethnic-specific susceptibility to disease.

A similar approach was used recently to identify genetic factors in animal models of diabetes and other complex human diseases, including Alzheimer’s disease[95]. A study by Zhong et al[68] (2010) in adipose tissue of C57BL/6-ob/ob × BTBR-ob/ob mice F2 progeny identified a strong causal subnetwork for T2D traits (called the “purple” module, enriched for genes involved in plasma glucose and insulin levels). They found that 37 eSNPs of genes in this module showed significant association with T2D in a GWAS report. Through additional prioritization steps and subsequent function validation studies, they identified mallic enzyme (ME1) as a key causal gene in this T2D subnetwork. A strong cis-eSNP of ME1 was associated with T2D. Future applications of such integrative genomic strategies in T2D or related disorders in human populations may prove insightful.

EQTL ANALYSES TO IDENTIFY GENE-ENVIRONMENT INTERACTIONS RELEVANT FOR T2D

As discussed above, GWAS have identified DNA sequence variants in the susceptibility to T2D, but these variants account for only a part of the estimated heritability[13,14]. Interactions between sequence variants and environmental stimuli are a logical step in better understanding the development of T2D. Thus, some of the missing heritability for T2D susceptibility may be explained by studies of the interaction between environmental factors and genetic variants or gene-environment (GXE) interactions[99]. Modeling GXE interactions in clinical or epidemiological settings is challenging and costly, due to few validated tools for assessing exposure (including dietary exposure), the need for large sample sizes, and the heterogeneity of exposures in populations[100-103]. Environmental factors usually influence insulin resistance and T2D risk over long periods of time; thus, accurate assessment of long-term exposure is needed to identify GXE interactions. A recent series of studies by Patel et al[104-106] utilized data resources from the National Health and Nutritional Examination Survey and integrated GWAS and environment-wide association studies to identify environmental factors, genetic factors, and GXE interactions involved in T2D susceptibility. However, they noted several significant limitations of such epidemiological approaches in adequately addressing influence of genetic variations on differences in environmental response in human populations.

Environmental factors, including diet and derived metabolites, can influence phenotypes by modulating gene expression in several ways. Variations in responses to environmental factors among individuals, and how these responses predispose to metabolic and other disorders, have been recognized[107]. Genetic variants modulate the environmental factor-mediated transcriptional response, which in turn dictates cellular response and may explain variability in metabolic responses to those factors[99]. Such dependency on external conditions or GXE interactions has been reported for genetic effects on gene expression in different organisms[108-110]. Transcripts responsive to environmental perturbation factors may manifest as eQTLs and are modulated by cis- and trans-eSNPs. A subset of these eSNPs associated with T2D, obesity, and/or glucose homeostasis traits may thus exhibit distinctive patterns of GXE eSNPs. Thus, identifying environmental factors that modulate insulin sensitivity and other early pathophysiological manifestations of T2D and its integration into eQTL analyses will further improve the power to construct causal gene regulatory networks involved in T2D susceptibility.

A few recent studies implemented a novel “cellular genomics” approach[111] to elucidate genetic controls on GXE interactions, critical to understanding the pathophysiology of complex diseases. In this novel paradigm, researchers analyzed the molecular consequences of genetic variants to assess interactions with environmental factors via quantification of processes (like gene expression) in cells from human subjects grown in uniform culture conditions. This concept is illustrated in Figure 5. Utilizing transformed lymphocytes, the studies examined genetic control in response to radiation, chemotherapeutic drugs, and hormones (glucocorticoids)[112-114]. Two similar studies in primary human cells mapped genetic regulators responding to growth factors (BMP-2), hormones (dexamethasone), cytokines (prostaglandin E2 in human osteoblasts), and oxidized low density lipoprotein (in human aortic endothelial cells)[115,116]. Despite the encouraging success of these studies, no studies so far have evaluated GXE interactions with a cellular genomic model relevant to T2D and related metabolic disorders. Although this model may miss some whole organism-level complexity[117] of T2D pathogenesis (which involves multiple tissues), it does represent an innovative approach by going from cellular to organismal phenotype analysis for identification of function of genetic variants involved in T2D susceptibility. Mapping GXE eSNPs for function-based prioritization of T2D and related metabolic disease-associated SNPs is a critical step towards designing efficacious strategies to reduce the public health burden of common metabolic disorders triggered by increased exposure to dietary and other environmental factors.

Figure 5 Types of gene-by-environment interactions in cellular genomic models to study gene-by-environment expression quantitative trait locis. Cells from a cohort of subjects are grown in pairs under uniform in vitro treated and untreated conditions to study environment-dependent or -independent effects of genotype on expression of transcripts (a quantitative trait). β1 and β2 are genotype effects on transcript expression under treated and control conditions, respectively. Different models of gene-by-environment (GXE) includes Null model: β1 = β2 = 0 (A and B); No-interaction eQTL model: β1 = β2 ≠ 0 (C); Treated-only expression quantitative trait loci model: β1 ≠ 0 and β2 = 0 (D); Control-only eQTL model: β1 = 0 and β2 ≠ 0 (E); and General interaction eQTL model: β1 ≠ 0 and β2 ≠ 0 but β1 ≠β2 (F). Black line indicates expression in cells under control condition (untreated) while blue line indicates expression in environmental/dietary factor treated cells.

EQTL AND PHARMACOGENOMIC STUDIES FOR T2D

Several classes of anti-diabetic medications are used for the treatment of T2D[118]. Pharmacogenomic studies reviewing the role of genetic variants on drug responses (including adverse drug reactions) have yielded significant findings, including novel disease mechanisms for several complex diseases[119]. But a similar success for T2D has not been achieved[5,120]. Pharmacological interventions using peroxisomal proliferator activated receptor gamma (PPARγ) agonists like pioglitazone improve insulin sensitivity and can reduce the risk of progression to T2D[121]. However, approximately 25% of patients do not respond adequately to PPARγ agonists[122]. Genome-wide transcriptomic analysis by our laboratory showed significant inter-individual variability in gene-expression response after pioglitazone treatment in people with impaired glucose tolerance[123]. However, little is known about the genetic architecture of variation in pioglitazone-mediated transcriptional response in human populations. Identifying the genetic variations that interact with pharmacological treatments like PPARγ agonists is of high clinical interest. eSNPs may modulate the expression of key transcripts in response to anti-diabetic drugs in target tissues and can explain the interindividual variability in treatment outcome[124,125]. Identifying genetic (and epigenetic) variants that modulate the pharmacological treatment-mediated transcriptional response, which in turn dictates the treatment outcome in T2D, is an open area of research. A novel approach that systematically characterizes the set of eSNPs involved in anti-diabetic medicine-mediated transcriptional modulation (gene-drug interaction eSNPs, or GXD eSNPs) in tissues relevant to glucose homeostasis will be useful in stratifying populations in efficacy studies, to improve the quality of clinical decision-making and treatment options for T2D.

FINDING EQTLS: END FOR A NEW BEGINNING

eQTL analyses provide statistical evidence for genotype-dependent variations in transcript abundance and should be considered a starting point for investigating the effects of DNA polymorphisms at the molecular level[34]. Transcript abundance depends on a dynamic relationship between transcript synthesis, stability, and degradation[48]. Thus, DNA polymorphisms may affect transcript abundance by several known and unknown mechanisms. Studies in human subjects have shown that sequence-specific regulation of mRNA expression is mediated by several molecular mechanisms, including allelic variability in transcription factor binding, chromatin remodeling, changes in DNase I hypersensitivity by histone methylation and acetylation, interaction between chromatin segments, alteration of splicing, sequence-dependent allele-specific DNA methylation, alteration of miRNA synthesis, and miRNA target binding[50-52,126-130]. GWAS-implicated variants for complex diseases are enriched in non-coding functional domains of the genome, including sequences involved in chromatin remodeling[131-133]. Many transcripts that show strong co-expression and cis-eSNPs for one transcript may appear as trans-eSNPs for a co-regulated transcript located in other chromosomes. Thus, a functional role of prioritized cis- and trans-eSNPs needs to be validated by appropriate molecular experiments to distinguish causal from correlative effects[134-136]. Studies have used allelic expression imbalance analysis, electrophoretic mobility shift assays, and transient transfection based luciferase reporter assays[56,137-141] to identify the molecular effects of genetic variants (cis-eSNPs) on gene expression; however, high-throughput methods are needed to validate in parallel the large number of findings from genomic studies[134,135,142]. Several novel high-throughput methods, including massively parallel reporter assays and massively parallel functional dissection, are now available to show evidence of causality for regulatory variants[143-146]. Functional relevance of the candidate eQTL transcripts in T2D pathophysiology also need to be validated by demonstrating their effects upon experimental up- or down-regulation in in vitro or in vivo experimental models[147,148].

In summary, many factors (including genetic, epigenetic and environmental factors) affect susceptibility to T2D. Instead of investigating different sources of variation in isolation, an integrative functional omics paradigm that traces early molecular changes through layers of biological information, including eQTLs, promises to be a useful approach[136]. Such an approach will promote optimal understanding of the etiology of T2D and lead to the identification of ethnic-specific primary causal variants. Ultimately, the knowledge gained from studies using these approaches can be used to build better classifiers of T2D risk than those based on DNA sequence variants alone.

ACKNOWLEDGMENTS

We thank Karen Klein (Translational Science Institute, Wake Forest University Health Sciences) and Ethel Kouba (Internal Medicine, Endocrinology) for critical reading and editing of our manuscript. We dedicate this article to the memory of Dr. Steven C. Elbein.