Transcriptomics in type 2 diabetes: Bridging the gap between genotype and phenotype.

Type 2 diabetes (T2D) is a common, multifactorial disease that is influenced by genetic and environmental factors and their interactions. However, common variants identified by genome wide association studies (GWAS) explain only about 10% of the total trait variance for T2D and less than 5% of the variance for obesity, indicating that a large proportion of heritability is still unexplained. The transcriptomic approach described here uses quantitative gene expression and disease-related physiological data (deep phenotyping) to measure the direct correlation between the expression of specific genes and physiological traits. Transcriptomic analysis bridges the gulf between GWAS and physiological studies. Recent GWAS studies have utilized very large population samples, numbering in the tens of thousands (or even hundreds of thousands) of individuals, yet establishing causal functional relationships between strongly associated genetic variants and disease remains elusive. In light of the findings described below, it is appropriate to consider how and why transcriptomic approaches in small samples might be capable of identifying complex disease-related genes which are not apparent using GWAS in large samples.

Introduction

The value of transcriptomic analysis in genetic studies

The transcriptome is the complete set of expression products transcribed from the genome in a specified tissue or population of cells. It is commonly used to refer to the expression of genome-wide protein coding genes and we follow that convention here, although our thesis applies equally to the expression of non-protein coding genes. In this review we discuss the use of whole genome gene expression approaches (referred to here as transcriptomics) in searching for the genes underlying complex genetic diseases. For brevity and clarity we focus on one such metabolic disease, type 2 diabetes (T2D), and its effect in Mexican Americans, as a typical example of a well-studied complex genetic disease. While we focus on this as a single complex genetic disease, it should be noted that T2D is strongly correlated with other complex disorders such as obesity, insulin resistance, the metabolic syndrome, dyslipidemia, cardiovascular disease, and a variety of diabetic complications. Thus, this discussion has wide relevance to other complex genetic diseases. It is our contention that transcriptomic approaches can provide a valuable bridge between the genotype and the phenotype. Recent studies have highlighted several issues in connecting genotypic approaches, such as (1) genome-wide association studies (GWAS), (2) genome-wide sequencing (GWS), and (3) expression quantitative trait locus (eQTL) analysis, with physiological disease phenotypes. We propose that the gap between these genotypic approaches and resulting phenotypes may be bridged or narrowed by the use of transcriptomic approaches. The reason for this is that transcriptomics occupy a unique position between the genomic or genotypic level and the physiologic or phenotypic level. We note that eQTL also provide a potential link between genotype and phenotype but eQTL typically rely on the connection between genotypes and gene expression and thus may be considered an enhanced form of GWAS. Gene expression as RNA provides the direct raw material for protein expression, which acts on and within cells and tissues, influencing so-called endo-phenotypes, and intermediate phenotypes and which eventually are measurable as physiological changes at the whole body level

In this review we also address potentially confounding factors such as the effect of the phenotype, including the disease phenotype, on transcription and on the myriad mechanisms by which proteins undergo post-translational modification. Nevertheless, these factors do not seriously affect the underlying rationale. A second issue is that transcriptomic analysis does not provide a direct causal link between specific nucleotides in the genome and downstream phenotypes such as disease-related clinical traits. This is also a limitation of GWAS and eQTL studies. However, transcriptomic analysis provides a link between RNA transcripts, which are copied directly from specific regions of the genome, and downstream measurable phenotypes. In this sense they provide a bridge between genotype and phenotype. Causal evidence requires careful follow-up experiments in model cells, tissues, organisms and eventually in humans. Transcriptomic experiments provide valuable initial correlational evidence to guide these follow-up studies.

Background

The diabetes pandemic and the burden of diabetes in Mexican Americans

Type 2 diabetes (T2D) is a major public health issue in the US, particularly in minority groups including Mexican Americans, the major sub-group of the US Hispanic population. Over the past 20 years, we have examined the genetics of T2D, obesity, metabolic syndrome and their related traits in low-income Mexican American families in San Antonio, Texas. T2D is a complex cardio-metabolic disorder characterized by both insulin resistance and pancreatic β-cell dysfunction [17] and is associated with significant morbidity from the micro- and macro-vascular complications of the disease. The compound burden of an increasing global T2D epidemic together with its co-morbid conditions including obesity, dyslipidemia and hypertension makes this a major world-wide public health problem [3]. The International Diabetes Federation estimates that globally 382 million people had diabetes in 2013, and that about 592 million people will be afflicted with diabetes by the year 2035. In the US, according to the Centers for Disease Control and Prevention [76], it is estimated that approximately 29.1 million people were afflicted with diabetes during the years 2010–2012, of whom about 21.1 million had diagnosed diabetes, while the remaining 8.1 million were unaware of their disease. The estimated prevalence of diagnosed diabetes based on age-adjusted data for individuals aged 20 years or older exhibits remarkable ethnic disparities (CDC, 2010–2012): non-Hispanic European Americans = 7.6%, African Americans = 13.2%, Asians = 9.0%, and Hispanics = 12.8% (of whom Mexican Americans = 13.9%).

Genetic variation in T2D and other complex diseases

T2D is a common, multifactorial disease that is influenced by genetic and environmental factors and their interactions. There are many large-scale ongoing efforts to localize and characterize T2D susceptibility genes using genome-wide association study (GWAS) approaches. To date, the GWAS method has achieved substantial success in localizing novel T2D susceptibility loci and loci for T2D-related glycemic traits (about 90 loci), obesity loci (~ 90), and loci for metabolic syndrome or its components (~ 50 loci), e.g. reviews: [4], [20], [28], [29], [41], [47], [51], [64], [65], [67]. However, common variants identified by GWAS explain only about 10% of the total trait variance for T2D [29] and less than 5% of the trait variance for obesity [4], indicating that a large proportion of trait variance, and of heritability, is still unexplained [48]. Therefore, there has been increased interest in the potential role of rare variants in the etiology of common complex diseases (rare variant/common disease hypothesis). These rare variants are postulated to have larger penetrance with potential functional consequences and could contribute to missing heritability [13], [26], [48]. Recent advances in next generation sequencing (NGS) technologies have made it possible to obtain complete information on rare and common variants across the whole exome, genome or within a targeted region [25] and new information on low frequency or rare variants, including exome chip based data, that are associated with T2D is emerging [16], [24], [45], [66], [73].

Functional translation of T2D association signals and the case for deep phenotyping

A large number of identified T2D loci has been associated specifically with defects in insulin secretion, highlighting the importance of pancreatic β-cell development and dysfunction. Beta cell function is commonly measured by the oral disposition index and the hyperglycemic clamp, which measures insulin secretory capacity. In contrast, relatively few loci have been discovered relating to insulin resistance which is commonly assayed by measuring levels of insulin and glucose during the oral glucose tolerance test (OGTT) and the gold standard euglycemic hyperinsulinemic clamp. However, recent studies in nondiabetic subjects have focused on loci associated with insulin resistance, e.g. reviews by [29], [47], [55], [68]. Since the overlap between the GWAS association signals with T2D and fasting glucose is partial, differential genetic influences on physiological and pathophysiological variation in glucose homeostasis are possible [29], [39]. These localization signals have, so far, yielded few causal gene identifications and the functional relevance of a majority of the implicated genetic variants has yet to be established. More than 85% of GWAS-identified variants, including those for T2D, fall in non-coding regions of the genome, highlighting their potential role in gene regulation [21], [22], [34], [39]. Efforts to functionally assess GWAS susceptibility loci are challenging. To enhance our knowledge of biological mechanisms and translational possibilities corresponding to specific genetic findings, the genotype-driven approach, aided with deep phenotyping, appears to be a more appealing and powerful strategy [52], [69]. Together with more precise characterization of clinical and pathophysiologic phenotypes, a range of deep physiological phenotypes (sometimes referred to as “endophenotypes” or intermediate phenotypes) and “omic-metrics” such as epigenomics, transcriptomics, proteomics and metabolomics can be examined in an integrated fashion to understand the molecular mechanisms underlying the phenotypic expression of diseases such as T2D [31], [40], [43], [63]. It is hoped that these more functionally-specific phenotypes will exhibit higher genetic signal-to-noise ratios and speed causal gene identification.

Role of transcriptomics in disease gene identification

Transcriptomic analysis, the measurement of gene expression at the genome-wide level, bridges the gulf between GWAS and physiological studies (Fig. 1). This is because of the causal relationships between the genome, the transcriptome and physiological traits. Information normally flows from the genome to traits of interest via the transcriptome. Thus, the information residing in genes predisposing to disease is often differentially expressed as RNA whose concentration can be correlated with disease-related traits at the physiological level. Recent GWAS studies have utilized very large population samples, numbering in the tens of thousands (or even hundreds of thousands) of individuals, yet establishing causal functional relationships between strongly associated genetic variants and disease remains elusive [4]. Finally, with the advent of inexpensive and accurate Whole Genome Sequencing (WGS) we are now on the cusp of a revolution in genomic association studies, involving genetic datasets of trillions of nucleotide measurements in combination with multiple physiological measurements [40]. As we enter this exciting era of genetic research, it may be an appropriate time to reassess some of the fundamental assumptions underlying these types of studies.

Fig. 1

Simplified depiction of the relationship between complex diseases, gene expression and the genome. Both GWAS and transcriptomics seek to deduce the underlying causal genetic relationships between DNA sequence variation and physiological variation. Transcript profiling analyzes a more direct and greatly simplified connection between the genome and complex disease pathology than GWAS. eQTL explicitly analyzes the relationship between DNA sequence variation and gene expression although their extension to phenotypic analysis can be performed in a similar manner to GWAS. It should be noted that the relationships in each case, between nucleotides and clinical traits, are correlative rather than causal. For clarity, the role of active non-coding RNAs is not included although their role is implicit. Their incompletely understood function would presumably mimic that of the relatively well-understood protein coding RNAs. Furthermore, their concentration can be accurately assayed by the same transcriptomic approaches used for protein-coding RNAs and their similar transcript variants.

GWAS vs physiological studies

In striking contrast to recent GWAS results is the remarkable paradox that many physiological correlative studies can be well-powered to detect clinically meaningful and scientifically generalizable information with modest sample sizes, often including less than 100 individuals, e.g. [11], [14], [15], [46], [54], [61]. The paradox lies in the fact that the physiological traits under investigation are often highly heritable and strongly genetically determined. The majority (> 85%) of GWAS signals are located in non-protein-coding DNA sequences, which collectively comprise 98.5% of the genome, and it is likely that their effect is due to their location within or near to regulatory regions which modulate gene expression. It is also likely that, given the GWAS data and what is currently known about the regulation of gene expression, multiple variants can act in concert to regulate the expression of genes which may be involved in key physiological processes, thereby contributing to complex multifactorial traits/diseases. Transcriptomics provides a method to detect these combinatorial effects by measuring changes in the “genetic phenotype” (i.e. gene expression) which is functionally intermediate between the genome sequence and clinically important traits. This biological gap is further reduced when deep physiological traits are used. By definition, these traits lie closer to gene effects than truncated clinical traits, such as the presence or absence of T2D or obesity, which are commonly used in large-scale genetic epidemiological investigations.

Complexity of gene regulation

Gene regulation is a complex multi-layered process involving numerous proteins and non-coding RNAs which may act at a great distance from their target gene. Elaborate multi-protein/RNA complexes must be assembled at the site of regulation. The regulatory mechanism may be intricate and variable, potentially involving transcript rearrangement and mRNA degradation. It is now clear that RNA has a diverse set of functions and is more than just a messenger between gene and protein. The mammalian genome is extensively transcribed, giving rise to thousands of RNA transcripts that are never translated into proteins. Whether all of these transcripts are functional is currently debatable, but it is evident that these include families of RNA molecules with a regulatory function [34]. The presence of a gene expression change, which is strongly correlated with relevant physiological changes, in the absence of proximate significant GWAS signals, suggests that relatively distant regulatory variants (and potentially many such variants) may act in combination to regulate the expression of the target gene of interest. Such putative gene expression-modulating variants could potentially act upon target gene expression through the mediation of non-protein-coding regulatory RNAs. For example, recent studies have shown that the expression of many genes is modulated by small interfering RNAs (siRNAs) and micro-RNAs, e.g. reviews by [10], [30], which do not encode proteins. In addition to microRNAs, many non-protein-coding RNA species (or “RNA genes”), such as long non-coding RNAs [42], are transcribed from the genome. Thus, there is compelling evidence that most of the genome may be transcribed [5], [6], [9], [19], [38], [53], [58], [59], [62] and the potential role of non-protein coding RNA genes in the modulation of protein-coding gene expression remains to be fully evaluated.

Mode of action of expression regulatory variants

The putative gene-expression-modulating variants could, in principle, act either through non-coding RNAs or through protein coding RNAs (via proteins). siRNAs have been shown to operate by binding to mRNAs, leading to premature mRNA degradation and hence down-regulation of gene expression [10]. Those proteins known to be involved in gene expression modulation operate by binding, in a coordinated sequential fashion, to enhanceosomes, repressosomes, spliceosomes and transcriptosomes (for review see [2]). In some cases the regulatory binding site has been shown to reside far from the putative target coding sequence, in regions between genes — so-called “gene deserts” [60], [72]. Regulatory binding sites thus integrate the effect of multiple signals and the effect of individual proteins or RNAs may be multiplicative, by analogy with the amplification effects of enzyme regulatory complexes found in intermediary metabolism and other cellular processes. Thus, such putative genetic mechanisms could evade GWAS approaches designed to detect the association of common or rare variants with whole-body physiological traits such as the presence or absence of T2D or obesity. In contrast, transcriptomic techniques measure gene action directly and can bypass the relatively weak signals of widely scattered variants within the genome.

Transcriptomics vs expression quantitative trait loci (eQTL) approaches

The transcriptomic approach uses quantitative gene expression data and disease-related physiological data to measure the direct correlation between the expression of specific genes and physiological traits. A related but conceptually distinct technique measures the correlation of genetic variants with gene expression. This approach seeks to identify genomic regions, termed expression quantitative trait loci (eQTL) that are significantly associated with the expression of specific genes. Thus, eQTL analysis lies between GWAS and the transcriptomic gene profiling described here. In some respects it may be considered a subtype of GWAS or an enhanced form of GWAS. However, unlike GWAS, eQTL analysis tends to require relatively smaller sample sizes. The power for a given sample size depends on the effect size and it is reasonable to expect larger genetic effects on transcription than on typical GWAS signals at the nucleotide level. Thus, sample size requirements may often be smaller than for GWAS. The number of statistical tests also affects power, since it determines the false positive error rate. In particular, the GWAS approach can be problematic because the massive number of statistical tests, between genetic variants and traits of interest, presents an enormous potential for false positive results. Whether considering genotyped variants of several hundred thousand or several million, including imputed variants, the scale of the problem is vast. In the case of whole genome sequencing, the problem is at least an order of magnitude greater. This also presents a problem for eQTL analyses, which measure the correlation between a similar number of genetic variants and the expression of transcripts derived from microarray or RNAseq data. In contrast, multiple testing presents much less of a problem for transcriptomic analysis, which measures the correlation between the expression of a comparatively small number of transcripts and physiological traits of interest.

Since eQTL signals involve a necessary intermediate step between genomic variants and physiological traits, they would be expected to provide support for some GWAS signals. This has been shown empirically to be the case with many loci for T2D [32]. The great majority of eQTL identified to date have been found to act in cis with the gene of interest rather than in trans [27]. Part of the explanation is the greatly reduced power to detect trans effects given the enormous multiple testing issues related to searching for trans eQTL, where most of the 3 billion bp of the genome must be searched, rather than a few hundred kbp for cis eQTL. The situation may be worse if there are multiple trans eQTL for a given gene of interest, scattered throughout the genome, and they act combinatorially and multiplicatively on gene expression (see Section 2.6 below). This may provide a partial explanation for the fact that many GWAS signals appear to act in trans. This is because a significant fraction of variants cannot, with certainty, be ascribed to act upon any particular gene because they tend to lie at some distance from any protein coding sequence and, therefore, their mode of action is less likely to occur in cis with respect to genes and potentially more likely to act in trans. A large proportion (43%) of GWAS hits are located in intergenic regions [33]. The lack of close proximity of many potential regulatory variants to coding sequences may in part be due to the fact that protein coding sequence accounts only for approximately 1.5% of the genome [35]. Nevertheless this assumption of major trans effects remains to be validated by rigorous laboratory testing of the effect of potential regulatory variants on the expression of hypothetical target genes. In general, such tests are not currently feasible on a large scale.

Gene expression variants and GWAS

In light of the findings described below (Section 2.7), it is appropriate to consider how and why transcriptomic approaches in small samples might be capable of identifying complex disease-related genes which are not apparent using GWAS in very large samples. With regard to the considerations discussed above for gene identification in complex diseases, we may consider a hypothetical simplified example. If two genetic variants exist which act combinatorially and multiplicatively to strongly regulate the expression of a protein-coding gene “G” but are located sufficiently far apart within the genome, such that the effect of one or both of the variants is operating “in trans” with respect to gene “G”, then the association of the two variants with a clinical trait of interest might be effectively invisible to detection by GWAS, even if both variants could be effectively targeted by specific oligonucleotides. This is because the association of each variant with the trait of interest may be too weak to reach significance in association testing even though the combined and multiplicative effects of the variants have a strong effect on the expression of “G”, which in turn influences the clinical trait. As previously noted, many GWAS signals are in fact located in non-exomic regions of the genome and are suspected of being regulatory variants or of tagging local regulatory sequences which are often non-canonical, such that a clear functional sequence pattern is difficult to discern. The situation becomes even more difficult if multiple genetic variants are involved (see below for a more complex hypothetical example) and if many or all of the variants are rare. In such a scenario it is not unreasonable to propose that association of gene expression signals with a clinical trait of interest may tend to be stronger than the association of individual nucleotide variants. As one example of such a putative mechanism, in a recent study of eQTL contributing to the risk for coronary artery disease, multiple cis and trans eQTL were identified, none of which contributed to the well-known 9p21.3 risk locus, containing the “first common variant associated with coronary artery disease to yield to the GWAS approach” [12]. Similarly, other studies have shown the presence of multiple trans eQTL which “are far from any GWAS association signals and thus cannot be identified from the GWAS alone” [32]. This concept finds additional support in the empirical example described below (Section 2.7).

Multiple distributed regulatory expression variants

Here we expand upon the theoretical example from the previous section to a more complex example. Instead of a single regulatory SNP acting in cis, relative to “G”, the regulatory region might instead, and more realistically, be hypothetically distributed over perhaps 10 genetic variants throughout the genome, all acting in trans or in combination with one or more cis variants. If all 10 variants are involved in the regulation of the “G” transcript, then the potential GWAS signal, that might otherwise have resided in a single SNP, could be diluted by several orders of magnitude. This effect would be worsened if truncated high-level phenotypes are used (e.g. the presence or absence of T2D or obesity) rather than quantitative deep physiological phenotypes (e.g. measurements from the oral glucose tolerance test [OGTT] and the euglycemic hyperinsulinemic clamp). In effect, the use of an intermediate “genetic phenotype”, gene expression, which is closer to both clinically relevant clinical phenotypes and to the immediate effect of gene action, further increases the power to detect a clinical effect of gene action. This is because the gap between the genomic measurement and the end-phenotype is greatly decreased and multiple genetic effects from individual nucleotides converge on a single gene expression phenotype. Thus, this “concentrating” effect suggests that the design of many GWAS experiments may tend to dilute the effects they seek to elucidate. Another factor enhancing transcriptomic gene identification compared with GWAS is that it isolates and effectively concentrates gene activity within a specific tissue and metabolic state, whereas GWAS necessarily targets the whole genome at all times and states.

Empirical example of gene discovery by transcriptomics

For illustrative purposes, an example of the transcriptomic approach described herein is provided by the recent identification of a gene, alcohol dehydrogenase 1B (ADH1B) with unexpectedly robust gene expression correlations with multiple deep phenotypes related to T2D. Detailed results are provided elsewhere [74]. RNA was extracted from abdominal subcutaneous adipose tissue biopsies from 75 unrelated non-diabetic Mexican American subjects and genome-wide gene expression measured on the Illumina HumanHT-12v4 Expression BeadChip platform, containing 47,324 oligonucleotide probes per array (derived from RefSeq release 38 and other sources), providing genome-wide transcriptional coverage of well-characterized genes, gene candidates, and splice variants. Microarray data were deposited in the Gene Expression Omnibus (GEO) database under ID GSE64567. ADH1B expression was inversely correlated with all obesity/insulin resistance (OB/IR) related variables at the tissue and whole body levels: waist circumference (P = 2.8 × 10− 9); BMI (P < 0.0001); fasting plasma insulin (P < 0.001); the Homeostasis Model Assessment of Insulin Resistance, HOMA-IR [50] (P < 0.01); Matsuda Insulin Sensitivity Index [49] (P < 0.01); β-cell function (Oral Disposition Index [71] across 120 min of the oral glucose tolerance test (ODI1–120; OGTT) (P = 8.6 × 10− 4), and pre-T2D (impaired glucose tolerance and fasting glucose during the OGTT) (P = 8.6 × 10− 5); and from the euglycemic hyperinsulinemic clamp technique [18]: total glucose disposal (TGD; P = 8.6 × 10− 10); hepatic glucose production (HGP, P = 8.6 × 10− 4); and free fatty acid IR index (FFA-IRI; P = 1.3 × 10− 4). The TGD rate mainly represents skeletal muscle IR; HGP indicates hepatic IR and FFA IR indicates adipose IR. As expected, ADH1B protein expression was decreased with high BMI (P < 0.05, N = 12). Transcription results for ADH1B were confirmed using quantitative PCR. In summary, decreased adipose ADH1B expression was associated with increased obesity, whole body IR, liver IR, skeletal muscle IR, adipose tissue IR, and declining β-cell function. While proof of causality remains to be established, the finding of strong consistent correlations across an array of phenotypes, including several deep physiological phenotypes, in a modest sample size, increases the prior probability that ADH1B is involved in T2D causation and helps prioritize functional testing. A cautionary note is that some of these physiological traits are, to some extent, correlated with each other and the apparent correlation between gene expression and these physiological traits could be partially explained by a reduced number of physiological variables.

In addition to the transcriptomic data obtained from 75 individuals, described above, the ADH1B chromosomal region (4q22) contains GWAS and eQTL data from approximately 69,000 and 1569 individuals, respectively. There were no significant GWAS signals found for association with T2D or related phenotypes (BMI, FPG, FPI) within at least 200 kb of ADH1B using GWAS SNP data (http://www.type2diabetesgenetics.org/accessed 4/30/2015). GWAS T2D and related trait data were obtained from a meta-analysis of 69,000 individuals of European ancestry and covered a total of 3.1 million SNPs. Similarly, no significant eQTL were found for ADH1B in subcutaneous adipose or any of more than 60 tissues examined (http://www.gtexportal.org/accessed 4/30/2015). The eQTL data were derived from 1569 RNAseq samples and contained 8930 significant genes. As previously noted, the lack of localized GWAS and eQTL signals should not be surprising and provides some support for the proposed model of gene action.

Relative gene expression and tissue availability for transcriptomic studies

The case for transcriptomic studies is likely to be strengthened when the gene of interest is highly expressed in the tissue(s) of interest and can, therefore, provide a potentially stronger and more accurate signal. According to the EMBL-EBI Expression Atlas (http://www-test.ebi.ac.uk) using data obtained by microarray assay of 76 normal human tissues (http://BioGPS.org; accessed 5/12/2015), ADH1B was most highly expressed in adipose tissue, with almost double the expression level found in the next highest tissue, liver. This finding was confirmed independently by high resolution sequencing of RNA (RNASeq) from 16 human tissues (Illumina Body Map 2.0; accessed 5/12/2015). In the RNAseq study, ADHIB had substantially higher expression in adipose tissue than in 15 other major tissues examined. Its expression, measured by RNAseq, was more than 3-fold higher in adipose tissue than in the next highest tissue, liver. On the other hand, transcriptomic studies in general are dependent on a limited source of available tissue which may not always be optimal for the studies' objective. In the present example, the tissue under study was highly relevant to the disease of interest and the results were not confounded by the presence of disease [27].

Some advantages and disadvantages of transcriptomic studies

A key advantage of transcriptomic studies is the ability to detect causative genes for complex genetic traits due to the proximity of transcripts and phenotypes of interest. A second potential advantage is the opportunity to measure environmental effects on transcription. Many transcriptomic studies tend to have case/control designs because of the difficulty of obtaining repeated samples of solid tissue, with the notable exception of liquid blood samples which can be easily sampled repeatedly over time. In either case, it is possible to incorporate environmental impacts in the study design, either by measuring differential environmental effects on the study groups or by measuring environmental effects over time and measuring the correlation of transcript expression with environmental exposure.

A limitation of transcriptomic studies is that relevant tissues are often not available and this difficulty is not easily overcome. One promising approach is the use of induced pluripotent stem cells (iPSCs) derived, for example, from human lymphoblastoid cell lines, although this too currently has limitations, including the accurate recapitulation of the target tissue phenotype. Secondly, some important genetic mechanisms may not be mediated through an effect on transcript levels, or protein levels, and thus may not be detectable by transcriptomic methods. Such mechanisms include nonsynonymous variants that may not change RNA or protein levels but may alter protein function. In this regard, transcriptomic studies can be complementary to genetic association or whole genome sequencing approaches.

Methods for dealing with confounding effects

Transcription is a dynamic process and it can be challenging to separate confounding effects on gene expression. One commonly used approach, used in the empirical example described above, in which disease was not a factor, is to use covariates such as age and sex in multiple regression analysis. Causal models have been developed to address issues of potential confounding. To take a hypothetical example, we may observe a correlation between a trait, such as BMI, and a gene's expression level. The obvious question to ask is whether this correlation, regardless of significance, is causal or not. If it is causal, then the gene's expression level will influence BMI, perhaps with some time lag. Conversely, there could be reverse causality, with the effect being reversed, such as the influence of BMI on a gene's expression level. More generally, there could be confounding of apparent causality by other unknown factors. As a hypothetical example, people with higher BMI might tend to come from socioeconomically deprived backgrounds that predispose them to eat low quality calorically dense foods more often, and this diet might also influence the expression of the gene(s) in question, even if the gene's expression itself does not lead to obesity

To disentangle such complex relationships, analytical methods have been developed, such as mediation analysis [44] and Mendelian randomization [8]. For example, the latter method works because genetic variants are randomized among us, by the rules of Mendelian inheritance (and some additional assumptions such as non-assortative mating) and are generally not correlated with possible confounders, such as someone's socioeconomic status and/or diet. Mendelian randomization is one of the valuable features that distinguish genetic epidemiology from most other fields of epidemiology. In the present case, if we have a genetic variant that is known to influence the expression level of the gene(s) in question, such as a cis-regulatory variant, then we could test whether this variant is also associated with BMI. If the gene's expression level causally influences BMI, then we should be able to observe an association of the cis-eQTL(s) with BMI as well. Alternatively, if causality is reversed, i.e. BMI influences gene expression, then there should be no such association between cis-eQTL(s) and BMI. This can be carried further by looking at effect sizes, e.g. if a doubling of gene expression correlates with a BMI increase of 10 units, then a cis-eQTL increasing gene expression by 50% should lead to a BMI increase of roughly 5 units. In any case, the use of this analytical strategy requires that we have cis-eQTLs for the gene(s) in question, or trans eQTLs, but this is much more unlikely. These eQTLs should come from the same tissue and preferably the same ethnicity.

Epigenomics, metabolomics and proteomics

In the context of gene expression modulation, three other related mechanisms should also be considered. (1) epigenomics/chromatin remodeling, (2) metabolomics, and (3) proteomics. Epigenomics and chromatin remodeling refer, respectively, to reversible heritable changes in DNA sequence and changes in chromosome structure. These occur most commonly via nucleotide and histone methylation and acetylation which affect localized chromosomal regions and gene expression via availability of DNA sequence and potential steric hindrance of transcriptosome formation. Epigenomics and associated chromosomal changes must be considered a potential mechanism of action in gene expression modulation. Epigenomic regulation of gene expression represents a sub-class of gene expression modulation and is, in principle, subject to the same kinds of mechanisms and processes as the modulation of RNA concentrations which are measured by transcriptomic profiling. However, as with GWAS's, whole genome epigenetic changes also occur at the level of the genome sequence and suffer from a similar difficulty in quantitation at the genome-wide level, and for similar reasons. Regions of epigenetic changes are influenced by multiple extrinsic and intrinsic factors and effectively represent changes at the genomic DNA level which can be considered analogous to DNA sequence differences similar to those measured by GWAS but with a greater variability, making epidemiological-scale measurements difficult to interpret [1], [70]

Complexity of proteomic and metabolomic analysis

Proteomics holds conceptually greater promise as an additional measure of disease-related genetic change given that it takes the “genetic phenotype” of RNA expression at least one step closer to the ultimate clinical phenotype of interest. However, with proteins and their multiple derivatives an enormous amount of complexity is introduced and this currently is difficult to resolve [75]. Proteomic and metabolomic measurements can be exquisitely sensitive to changes in the concentrations of specific proteins and metabolites [23]. However, the task of de novo identification of novel protein products and metabolites is more daunting than that of RNA profiling. In comparison, transcriptomic measurements, which rely on the identification of simple four-nucleotide polymers, whether measured by oligonucleotide arrays or RNA sequencing, are currently more amenable to disease gene identification. Thus, genome-wide analyses of very large numbers of proteins and metabolites, potentially orders of magnitude greater than RNA species, currently are limited by their inability to simultaneously identify and measure all possible species encoded by the genome. This is a technical issue which eventually should be resolved by better high-throughput high-resolution technology and likely will become a valuable adjunct to future transcriptomic analysis

It should be noted that there is not a simple 1:1 relationship between RNA expression and protein expression (including expression of covalently modified proteins) and this situation is even more complex for the biosynthesis and use of metabolites. This lack of a simple relationship between RNA expression and protein/metabolite expression makes interpretation of proteomic and metabolomic data challenging. In this regard, one lesson that may be drawn from GWAS analyses is that larger datasets may not always provide simple answers in multifactorial diseases where the underlying structure is very complex. For example, a small proportion of a given protein may become highly active after phosphorylation and become capable of initiating an exponential intracellular catalytic cascade of activity out of proportion to the expression of the original non-phosphorylated protein. Thus, the concentrations of multiple phosphorylated and dephosphorylated proteins may need to be measured (not to mention the concentrations of numerous other types of covalently modified proteins). As with GWAS, eQTL and transcriptomic correlations, the measurable proteomic and metabolomic changes are ultimately encoded within the genome and transferred via expression of RNAs to eventually interact with the intracellular and extra-cellular environment.

Identification of accessible chromosomal regions for gene expression regulation

Newer epigenomics-related approaches can aid the search for meaningful non-coding variation by focusing on regions of the genome where regulatory protein complexes are likely to be bound, so-called nuclease hypersensitivity sites. Such methods are based on the location of genome-wide sterically accessible open chromatin regions. These techniques exploit the power of massively parallel next generation sequencing (NGS) coupled with nuclease digestion to elucidate the genome-wide location and sequences of gene expression regulatory sites. Current methods include: (1) the assay for transposase-accessible chromatin using sequencing (ATAC-seq) [7], which locates putative sites for all DNA-binding proteins, and (2) chromatin immunoprecipitation with sequencing (ChIP-seq) [36], [37] which locates the binding sites of known or suspected DNA-binding proteins. We have previously used ChIP-seq to elucidate genome-wide putative binding sites for the TCF7L2 transcription factor protein and identified binding sites for TCF7L2 regulation of expression of hepatocyte gluconeogenic genes which appear to be involved in hepatic insulin resistance [56], [57]. Although ATAC-seq and ChIP-seq analyses can be performed independently of transcriptomic analysis, they represent powerful synergistic approaches to identify specific sites of regulation of gene expression uncovered by transcriptomic studies.

Conclusion

In this review we have presented evidence that transcriptomic studies may offer insights into physiological mechanisms that are not easily obtained from GWAS or eQTL. We propose that the effects of multiple genetic variants may “converge” on a single genetic mechanism that can be assayed by calculating the correlation between the levels of a particular transcript and phenotype, whereas the effects of individual genetic variants are comparatively small and difficult to detect. This approach could allow transcriptomic studies to yield insights from much smaller sample sizes than have hitherto proven to be necessary for GWAS, particularly if “deep” phenotypes, that are closer to the physiological traits of interest, are analyzed. We discussed an empirical example where highly significant associations were detected between the levels of the transcript for alcohol dehydrogenase 1B, measured in adipose tissue, and multiple measures of obesity and insulin resistance. We conclude that transcriptomic analysis can provide a more convenient and direct path bridging the genotype and phenotype than is currently available using either GWAS or eQTL approaches.

Conflict of interest statement

The authors warrant that they have no conflict of interest with regard to the material in this manuscript.

Author agreement

All authors have seen and approved the manuscript being submitted. They warrant that the article is the authors' original work, has not received prior publication and is not under consideration for publication elsewhere.