Genetic effects on gene expression across human tissues.

Characterization of the molecular function of the human genome and its variation across individuals is essential for identifying the cellular mechanisms that underlie human genetic traits and diseases. The Genotype-Tissue Expression (GTEx) project aims to characterize variation in gene expression levels across individuals and diverse tissues of the human body, many of which are not easily accessible. Here we describe genetic effects on gene expression levels across 44 human tissues. We find that local genetic variation affects gene expression levels for the majority of genes, and we further identify inter-chromosomal genetic effects for 93 genes and 112 loci. On the basis of the identified genetic effects, we characterize patterns of tissue specificity, compare local and distal effects, and evaluate the functional properties of the genetic effects. We also demonstrate that multi-tissue, multi-individual data can be used to identify genes and pathways affected by human disease-associated variation, enabling a mechanistic interpretation of gene regulation and the genetic basis of disease.

The human genome encodes instructions for the regulation of gene expression, which varies both across cell types and across individuals. Recent large-scale studies have characterized the regulatory function of the genome across a diverse array of cell types, each from a small number of samples[1-3]. Measuring how gene regulation and expression vary across individuals has further expanded our understanding of the functions of healthy tissues and the molecular origins of complex traits and diseases[4-9]. However, these studies have been conducted in limited, accessible cell types, thus restricting the utility of these studies in informing regulatory biology and human health.

The Genotype-Tissue Expression (GTEx) project was established to characterize human transcriptomes within and across individuals for a wide variety of primary tissues and cell types. Here, we report on a major expansion of the GTEx project that includes publicly available genotype, gene expression, histological and clinical data for 449 human donors across 44 (42 distinct) tissues. This enables the study of tissue-specific gene expression and the identification of genetic associations with gene expression levels (expression quantitative trait loci, or eQTLs) across many tissues, including both local (cis-eQTLs) and distal (trans-eQTLs) effects.

In this study, we associate genetic variants with gene expression levels from the GTEx v6p release. We found pervasive cis-eQTLs, which affect the majority of human genes. In addition, we identify trans-eQTLs across 18 tissues and highlight their increased tissue specificity relative to cis-eQTLs. We evaluate both cis- and trans-eQTLs with respect to their functional characteristics, genomic context, and relationship to disease-associated variation.

Study design

The GTEx project has created a reference resource of gene expression levels from ‘normal’, non-diseased tissues. Every tissue sample was examined histologically; the sample was accepted for the project if the tissue was non-diseased and in the normal range for the age of the donor. RNA was isolated from postmortem samples in an ongoing manner as donors were enrolled into the study. For this data release, 44 sampled regions or cell lines were considered, each from at least 70 donors, and thereby considered suitable for eQTL analysis: 31 solid- organ tissues, 10 brain subregions including duplicates of two regions (cortex and cerebellum), whole blood, and two cell lines derived from donor blood and skin samples. We hereafter refer to these tissues, regions, and cell lines as the ‘tissues’ used in eQTL analysis. A total of 7,051 samples from 449 donors represent the GTEx v6p analysis freeze (Fig. 1a; Supplementary Information 1–5; Supplementary Figs 1–6; Supplementary Tables 1–10). RNA sequencing (RNA-seq) samples were sequenced to a median depth of 78 million reads. This is 4.3 times more samples than reported in the GTEx pilot phase[10]. DNA was genotyped at 2.2 million sites and imputed to 12.5 million sites (11.5 million autosomal and 1 million X chromosome sites) using the multi-ethnic reference panel from 1000 Genomes Project Phase 1 v3[11]. Sampled donors were 83.7% European American and 15.1% African American. Whole-genome sequencing was performed for 148 donors to a mean coverage greater than 30×, and all donors were exome-sequenced to a mean coverage over captured exons of 80×. The resulting data provide the deepest survey of individual- and tissue- specific gene expression to date, enabling a comprehensive view of the impact of genetic variation on gene expression levels. All data are available from dbGaP (accession phs000424.v6.p1) with multiple data views publicly available from the GTEx Portal (www.gtexportal.org).

Figure 1

Sample size and eGene discovery in the GTEx v6p study

a, Illustration of the 44 tissues and cell lines included in the GTEx v6p project with the associated number of cis- (left) and trans-eGenes (right) and sample sizes. Each tissue has a unique colour code (defined in Supplementary Fig. 5). b, Fraction of genes that are eGenes across all tissues by transcript class. The three tissues highlighted are: testis, which has the highest proportion of trans-eGenes; skeletal muscle, which has the largest sample size; and fibroblasts, which have the highest proportion of cis-eGenes. Dark bars depict the fraction of all curated human genes in GENCODE v19. Light bars depict the fraction of genes expressed in one or more tissues. c, Proportion of expressed genes that are cis-eGenes (y-axis) as a function of tissue sample size (x-axis). Colours represent tissues, as in a. d, Number of trans-eQTLs (x-axis) per tissue (y-axis), with sample size indicated by point size.

Expression QTLs across human tissues

We identified associations between the expression levels of all expressed genes (eGenes) and genetic variants (eVariants) located within 1 Mb of the target gene’s transcription start site (TSS), which we refer to as cis-eQTLs for convenience, without requiring evidence of allelic effects at each locus. However, the majority of cis-eQTLs do exhibit allele specific expression. We applied a linear model controlling for ancestry, sex, genotyping platform and latent factors[12] in the expression data for each tissue that may reflect batch or other technical variables (see Methods; Extended Data Fig. 1, Supplementary Information 6 and Supplementary Figs 7–10). Considering all tissues, we found a total of 152,869 cis-eQTLs for 19,725 genes, representing 50.3% and 86.1% of all known autosomal long intergenic noncoding RNA (lincRNA) and protein-coding genes, respectively (Fig. 1a, b). We identified a median of 2,816 autosomal protein-coding or lincRNA eGenes at a 5% false discovery rate (FDR) within each tissue (Extended Data Fig. 2a). Protein-coding genes without a cis-eQTL in any tissue were more likely to be expressed at low levels or loss-of-function intolerant and were enriched for functions related to development and environmental response, indicating specific contexts in which additional eQTLs may be identified (Extended Data Fig. 3). To identify cis-eGenes affected by more than one functional regulatory variant, we applied forward–backward stepwise regression (see Methods). This approach identified an additional 24,886 secondary cis-eQTLs, with 41.2% of protein-coding genes and 24.8% of lincRNAs having multiple, conditionally independent eVariants in at least one tissue (Supplementary Fig. 11).

Extended Data Figure 1

Association of known covariates with expression components removed by PEER

Each cell depicts the adjusted (R2) between a pair of variables. Scale bar specific to each panel is displayed at the bottom. Grey cells represent pairs of variables without sufficient data to estimate correlation. a, For each tissue, adjusted (R2) reflecting the proportion of variance explained by each covariate of the entire PEER component removed from the expression data. A selected set of the most relevant sample-specific covariates is shown here. b, For each tissue, adjusted (R2) reflecting the proportion of variance explained by each covariate of the entire PEER component removed from the expression data. A selected set of the most relevant donor-specific covariates is shown here. See Supplementary Information for complete set of covariates. c, Adjusted R2 between each PEER factor and known sample covariates in skeletal muscle. d, Adjusted R2 between each PEER factor and known donor covariates in skeletal muscle.

Extended Data Figure 2

Controlling the discovery of cis- and trans-eGenes

a, b, The number of unique cis- (a) and trans-eGenes (b) across all tissues at varying FDR thresholds. c, The per cent change in the number of cis-eGene discoveries comparing FDR (Storey’s q value) calculated across all tested genes in each tissue to FDR calculated only across autosomal lincRNA and protein-coding genes. Results are shown for each tissue and are stratified by gene type. The per cent change for each tissue and gene type is calculated as 100 ×(no. of cis-eGenes from q value on the restricted gene set – no. of cis-eGenes from q value on all tested genes)/(no. of cis-eGenes from q value on all tested genes). d, The per cent change in the number of cis-eGene discoveries comparing q value and the Benjamini–Hochberg (BH) procedure applied to all tested genes in each tissue. The percent change for each tissue and gene type is calculated as 100 ×(no. of cis-eGenes from BH – no. of cis-eGenes from q value)/(no. of cis-eGenes from q value). Box plots depict the IQR, whiskers depict 1.5× IQR.

Extended Data Figure 3

Properties of cis-eGenes and non-eGenes

In these analyses, ‘non-eGenes’ refers to the set of genes with no significantly associated cis-eQTLs. a, Density of expression (mean across samples, median across tissues, in log10 scale) for cis-eGenes and non-eGenes. The difference in mean expression for cis-eGenes and non-eGenes is significant (Wald test; P <1 ×10−16). b, Box plots comparing the probability of being loss-of-function-intolerant (pLI) scores[41] for cis-eGenes and non-eGenes, stratified by expression percentile across all genes (mean across samples, median across tissues, in log10 scale). On average, highly expressed non-eGenes have higher pLI scores than highly expressed cis-eGenes (t-test for the difference in mean pLI score between cis-eGenes and non-eGenes; BH adjusted P = 8.3 ×10−4 and 2.1 ×10−9 for (50,75]% and (75,100]%, respectively) and lowly expressed non-eGenes (t-test for the difference in mean pLI score between highly and lowly expressed non-eGenes; BH adjusted P ≤7.8 ×10−46, P ≤5.5 ×10−17 and P ≤2.1 ×10−3 for (75,100]% vs [0,25]%, (75,100 vs (25,50]% and (75,100]% vs (50,75]%, respectively). c, Gene Ontology (GO) analysis of tested protein-coding non-eGenes. We used the PANTHER overrepresentation test[78] (release 20160715) against 20,972 human genes as background to test for enrichment in GO biological processes using the GO database release 2017-02-28. Significant GO IDs (Bonferroni adjusted P < 0.05) were selected for analysis with REVIGO[79] to group similar ontological terms, which yielded 22 over-represented GO IDs. d, GO analysis of a more stringent set of protein-coding non-eGenes. Selected genes included those not tested in GTEx (532 genes) or those with a minimum nominal P value across tissues greater than 0.1 (692 genes). Of these stringent 1,224 non-eGenes, 808 were mapped in the GO analysis. Using a similar approach as in c, 20 over-represented GO IDs were identified. For both c and d, the x-axis represents the −log10 P value resulting from GO analysis. GO IDs are coloured by the broader enrichment category to which each corresponds. Box plots depict the IQR, whiskers depict 1.5× IQR.

To identify trans-eQTLs, we tested for association between every protein-coding or lincRNA gene and all autosomal variants where the gene and variant were on different chromosomes. To minimize false positives in trans-eQTL detection, we controlled for the same observed and inferred confounders as in the cis-eQTL analysis, and further removed genes with poor mappability, variants in repetitive regions, and trans-eQTLs between pairs of genomic loci with evidence of RNA-seq read cross-mapping due to sequence similarity. Applying this approach, we identified 673 trans-eQTLs at a 10% genome-wide FDR. This includes 112 distinct loci (R2 ≤ 0.2) and 93 unique genes (94 total gene associations, including a trans-eGene detected in both testis and thyroid) in 16 tissues (Extended Data Table 1, Extended Data Fig. 2b, Supplementary Information 7, 8, Supplementary Figs 12–15 and Supplementary Table 11). An alternative approach to quantifying FDR at the gene level (Supplementary Information 8 and Supplementary Table 12) identified 46 genes at 10% FDR, with estimated q values of less than 0.4 for all 94 gene associations identified using the genome-wide FDR[16]. By investigating long-range intra-chromosomal eQTLs (≥5 Mb from the TSS), we discovered an additional 33 eGenes (10% FDR; Extended Data Table 2 and Supplementary Information 9). We found decaying support for cis-regulation (or interaction between cis- and trans-effects) over increased genomic distances based on evidence of allelic effects (Extended Data Fig. 4). Evidence of cis-regulation fell below background levels between 0.85 and 1.3 Mb from the TSS, empirically supporting the conventional distance threshold of 1 Mb for cis-eQTL detection.

Extended Data Table 1

trans-eVariant and trans-eGene discoveries for genome-wide FDr control, and trans-eGene discoveries for gene-level FDr control

		Genome wide		Gene-level FDR
Tissue	No. of samples	No. of trans-eGenes	No. of trans-eVariants	No. of trans-eGenes
Muscle – Skeletal	361	9	43	4
Whole Blood	338	1	2	1
Skin – Sun Exposed (Lower leg)	302	6	16	3
Adipose – Subcutaneous	298	2	7	0
Lung	278	2	2	2
Thyroid	278	21	181	3
Cells – Transformed fibroblasts	272	1	10	1
Nerve – Tibial	256	0	0	1
Esophagus – Mucosa	241	3	11	3
Artery – Aorta	197	1	1	1
Skin – Not Sun Exposed (Suprapubic)	196	1	1	2
Stomach	170	0	0	2
Colon – Transverse	169	2	10	2
Testis	157	35	267	16
Pancreas	149	2	12	1
Adrenal Gland	126	1	1	1
Brain – Putamen (Basal ganglia)	82	3	11	2
Vagina	79	4	27	1
Total unique		93	602	46

Each tissue with non-zero values is included as a row; the columns include the number of samples for that tissue, followed by the number of unique trans-eGenes and trans-eVariants identified in the genome-wide tests, and the number of unique trans-eGenes found using gene-level FDR calibration (Supplementary Information 8). Ultimately, the set of 673 trans-eQTLs identified in the genome-wide approach yielded 602 unique trans-eVariants.

Extended Data Table 2

Distal (>5 Mb) intra-chromosomal eQTLs across the GTEx tissues

Tissue	No. of samples	No. of trans-eGenes	No. of trans-eVariants
Whole Blood	338	4	17
Skin – Sun Exposed (Lower leg)	302	2	12
Adipose – Subcutaneous	297	1	2
Lung	278	2	37
Thyroid	278	5	17
Cells – Transformed fibroblasts	272	4	17
Esophagus – Mucosa	241	15	184
Testis	157	4	30
Brain – Putamen (Basal ganglia)	82	1	2
Brain – Hypothalamus	81	1	2
Total unique	7047	33	284

Only tissues with at least one distal intra-chromosomal eQTL are listed.

Extended Data Figure 4

Identification of cis-acting eQTLs using allele-specific expression at chromosome-wide distances

a, A logistic regression based model was developed to predict the probability of phasing error as a function of distance and variant minor allele frequencies. When applied to chromosome 2 of 1000 Genomes sample NA12878, this model had a receiver operating characteristic (ROC) area under the curve (AUC) of 0.87 using population phasing compared to transmission phasing. b, ROC when applying the beta-binomial mixture model to detect cis-acting regulation to the GTEx v6p subcutaneous adipose cis-eQTLs, with an AUC of 0.88. As the null, eGenes were shuffled with respect to eVariants. c, Power analysis using all nominally significant (P <1.0 ×10−5) linkage disequilibrium pruned associations within 100 kb of the TSS illustrating the number of eQTLs with nominally significant (P ≤ 0.01) evidence of cis-regulation as a function of phasing error and eQTL effect size. Expression QTL effect size was calculated using a companion method[28], and uniform phasing error between 0 and 100% was introduced in silico. d, Proportion of nominally significant (P <1.0 ×10−5) linkage disequilibrium-pruned intrachromosomal eQTLs with nominally significant (P ≤ 0.01) ASE supported evidence of cis regulation in bins of increasing TSS distance. Observed indicates what is seen in the data, while Max Error indicates what would be expected in the worst-case scenario of phasing error (50%). e, Example of significant ASE supported cis-regulation at a distance of 52.7 Mb between eVariant rs17494053 and eGene ENSG00000108509 in whole blood. Each point represents allelic imbalance in a single eVariant homozygote (circle) or heterozygote (triangle).

As expected, sample size greatly affects eQTL mapping. Discovery of eGenes increased with sample size with no evidence of saturation at the full sample size for each tissue (Fig. 1c). The tissue with the highest number of identified cis-eGenes was tibial nerve, with 8,087 eGenes in 256 samples. Testis had the most trans-eGenes, with 35 eGenes in 157 samples (Fig. 1d), consistent with the elevated number of expressed genes (16,853 protein-coding genes and 4,362 lincRNA genes) and cis-eGenes (6,796 genes). Continued discovery of eGenes with increasing sample size suggests that the expression of nearly all genes is influenced by genetic variation (Extended Data Fig. 5a, b). We further observed that, for sub-sampled data ranging from 70 to 250 donors, sample size was a more significant contributor than additional tissues to the discovery of novel cis-eGenes (Extended Data Fig. 5c). For trans-eQTL mapping, we used informed subsets of variants to reduce the number of tests by one to three orders of magnitude (Supplementary Information 10 and Supplementary Table 13). We found statistical power to detect additional associations in these restricted tests, such as the test restricted to cis-eVariants. Our results indicate that ongoing increases in sample size will continue to yield additional eQTLs, both in the cis-eQTL setting, where smaller and conditionally independent effects will be identified, and in the trans-eQTL setting, where statistical power is the main limitation.

Extended Data Figure 5

Effects of sample size and assayed tissues on cis-eGene discovery

a, Number of significant cis-eGenes at 5% FDR (y-axis) discovered in 44 GTEx tissues versus sample size (x-axis). b, Number of significant eGenes at FDR 5% (y-axis) discovered in nine GTEx tissues each subsampled to various sizes where possible (n =70, 100, 125, 150, 175, 200, 225, 250, 275, 300, 325, and 350; x-axis). We computed the number of cis-eGenes at each subsample size (circles connected by lines). We also plotted the number of cis-eGenes discovered with no subsampling of donors (diamonds). c, Number of significant cis-eGenes at FDR 5% (y-axis) as a function of sample size and number of tissues assayed (x-axis). Each tissue was subsampled to 70, 100, 150, 200, and 250 donors, and a forward search was used to assess sequential combinations of tissues that maximize the total number of unique cis-eGenes discovered.

Allele-specific expression across human tissues

The effect of cis-regulatory variation can also be quantified by allele- specific expression (ASE) analyses obtained by measuring the allelic imbalance of RNA-seq reads at transcribed heterozygous sites. A large-scale, multi-tissue resource of ASE estimates complements eQTL mapping by providing access to individual-specific effects, which assists in the interpretation of rare variants, somatic mutations and patterns of imprinting[8,13,14]. We measured ASE at more than 135 million sites across tissues and donors, with a median of over 10,000 genes quantified per donor (Supplementary Information 11 and Supplementary Fig. 16). In total, 63% of all protein-coding genes could be tested for ASE in at least one donor and tissue, with 54% exhibiting significant allelic imbalance (binomial test, 5% FDR, |effect size| ≥ 1, Extended Data Fig. 6a). Overall, 88% of testable genes had significant allelic imbalance in at least one donor (binomial test, 5% FDR), demonstrating an abundance of cis-linked regulatory effects. Per donor, a median of 1,963 genes had significant allelic imbalance in at least one tissue, with a median of 570 genes where the donor was not heterozygous for a top eVariant, suggesting more complex or rarer regulatory effects at these loci.

Extended Data Figure 6

Measuring cis-regulatory variation using ASE

a, Proportion of protein-coding genes with ASE data in at least one tissue as a function of donors (top row) and with significant imbalance (binomial test versus 0.5, 5% FDR) stratified by ASE effect size (|log2(hap count/ hap count)|) deciles. Gene-level measurements of haplotype expression were calculated by aggregating counts per sample across all heterozygous variants with ASE data within the gene using population phasing. The following filters were applied on ASE data: total coverage ≥8 reads, no mapping bias in simulations[76], UCSC mappability > 50, and no significant (FDR > 1%) evidence that variant is monoallelic in expression data[75]. b, log2 transformed cis-eQTL effect size (x-axis) versus log2 transformed ASE effect size (y-axis) for whole blood (Spearman’s ρ =0.82) and c, subcutaneous adipose (Spearman’s ρ =0.74).

Tissue-sharing and specificity of eQTLs

The extensive and diverse tissue sampling allowed us to develop a global view of how genetic effects vary between tissues of the human body by evaluating the sharing of eQTLs across tissues. We performed a meta-analysis across all 44 tissues for both cis- and trans-eQTLs to assess eQTL sharing between tissues. To do so, we applied Meta-Tissue[15], a linear mixed model that allows for heterogeneity in effect sizes across tissues and controls for correlated expression measurements that result from collecting multiple tissues from the same donors. For each eQTL, we estimated the posterior probability that the effect is shared in each tissue (m value). For both cis- and trans-eQTLs, we observed patterns that reflected relationships between related tissues and concordance between cis and trans in estimates of tissue similarity (Fig. 2a, Supplementary Information 12 and Supplementary Fig. 17). The strongest broad pattern observed was the high correlation among brain tissues (median Spearman’s ρ of 0.584 (cis) and 0.241 (trans)) and among non-brain tissues (median Spearman’s ρ of 0.606 (cis) and 0.165 (trans)), with much lower correlation observed between these two groups (median Spearman’s ρ of 0.499 (cis) and 0.096 (trans)). Within non-brain tissues, we observed strong correlation among closely related tissues, such as arterial tissues (median Spearman’s ρ of 0.743 (cis) and 0.264 (trans)), skeletal muscle and heart tissues (median Spearman’s ρ of 0.672 (cis) and 0.184 (trans)), and skin tissues (Spearman’s ρ of 0.804 (cis) and 0.365 (trans)). Overall, the median pairwise correlation between tissues was 0.547 (cis) and 0.138 (trans).

Figure 2

Patterns of tissue sharing of eQTL effects

a, Similarity (Spearman’s ρ) of Meta-Tissue effect sizes between tissues for cis- (upper triangle, 5% FDR) and trans- (lower triangle, 50% FDR) eQTLs. Tissues (by colours as in Fig. 1a) are ordered by agglomerative hierarchical clustering of the cis-eQTL results. b, The number of tissues in which a given eQTL is shared as a function of tissue sample size. For each tissue, we estimated the degree of sharing (number of tissues with m >0.9) for all eQTLs identified in that tissue at a 5% FDR. Tissues were then binned into quartiles on the basis of sample size. A higher proportion of eQTLs identified in tissues with small sample sizes have shared effects across multiple tissues compared with more deeply sampled tissues. This pattern inverts at higher sample sizes where more of the effects are tissue-specific. The median number of shared tissues is plotted for each quartile as a horizontal black line. c, Distribution of the number of tissues having Meta-Tissue m > 0.5 for the top variant for each trans-eGene at 50% FDR, and FDR-matched, randomly selected cis-eGenes (also 50% FDR). cis-eGenes were matched for discovery tissue to the trans-eGenes.

The patterns of sharing were also supported by replication between single-tissue cis-eQTLs, estimated by π1 (the proportion of true positives[16]) among the eQTLs identified in one tissue and then tested for replication in a second tissue (Extended Data Fig. 7a, median π1 = 0.740). The patterns held even when accounting for variable number of overlapping donors among pairs of tissues in the GTEx study design (Extended Data Fig. 7b–e). cis-eQTLs exhibited a distinctly bimodal pattern of tissue sharing—they were likely to be either shared across most of the 44 tissues or specific to a small subset of tissues (Fig. 2b). This bimodality was further supported by three different methods: simple overlap of the single tissue eQTLs, a hierarchical multiple comparison procedure (treeQTL[17]), and an empirical Bayes model (MT-eQTL[18]; Extended Data Fig. 8a–c). Each method also demonstrated that cis-eQTLs discovered in tissues with larger sample sizes were more often tissue-specific; however, estimates of tissue-specificity for large sample-size tissues can be influenced by difficulty in replicating small effect-size eQTLs in tissues with fewer samples.

Extended Data Figure 7

Replication of cis-eQTLs between tissues

a, π1 statistics for cis-eQTLs are reported for all pairwise combinations of discovery (y-axis) and replication (x-axis) tissues. Higher π1 values indicate a stronger replication signal. Tissues are grouped using hierarchical clustering on rows and columns separately with a distance metric of 1 − ρ, where ρ is the Spearman’s correlation coefficient of π1 values. π1 is calculated only when the gene is expressed and testable in the replication tissue. b, c, π1 replication is reported between tissues subsampled down to 70 non-matched (b) and matched (c) donors. In (c), grey tiles indicate tissue pairs with fewer than 70 shared donors. d, Effect of sample size (x-axis) on average π1 replication (y-axis) across all replication tissues. e, Scatter plot of π1 scores among tissue pairs with matched (x-axis) and non-matched (y-axis) donors.

Extended Data Figure 8

Tissue-specificity of cis- and trans-eQTLs

a, Sharing of independently identified cis-eGenes across the 44 GTEx tissues (cis-eGenes are independently identified in each of the 44 tissues and then binned by the number of tissues in which they appear). b, Sharing of cis-eGenes across 44 GTEx tissues that were identified using the hierarchical multi-tissue analysis. c, The prior probabilities of having significant variant–gene association in different numbers of tissues, calculated using an empirical Bayes model. The prior probabilities are summed up on the basis of the Hamming weights of the corresponding cis-eQTL configurations. d–g, Meta-analysis performed using Meta-Tissue for trans-eGenes (50% FDR), randomly selected cis-eGenes (50 % FDR), and an equal number of the top cis-eGenes by P value. Distribution of the number of tissues that have Meta-Tissue m values greater than a given threshold (d, 0.5; e, 0.6; f, 0.9) across variant–gene pairs that have an effect (based on m value thresholding) in at least one tissue. g, The same distribution as d except that variant–gene pairs with predicted effect in zero tissues (based on the number of m values > 0.5) are included. Meta-Tissue predicts that many cis-eGenes (50% FDR) and trans-eGenes (50% FDR) will have an effect in zero tissues. The number of zero predictions is largely reduced for the top most significant cis-eGenes. h, Distribution of observed replication rate between pairs of tissues for trans-eQTLs (10% FDR) versus the predicted replication rate for trans-eQTLs (10% FDR) based on a negative binomial generalized linear model trained on cis-eQTLs (10% FDR0.1). This model directly accounts for effect size and minor allele frequency. Replication rates shown for a range of FDR thresholds in replication tissue. Box plots depict the IQR, whiskers depict 1.5× IQR.

Overall, we observed much greater tissue specificity for trans-eQTLs than a set of FDR-matched cis-eQTLs (Fig. 2c); this observation was robust to choices of m value threshold and selection criteria for matching cis-eQTLs (Extended Data Fig. 8d–g). While 3.8% of trans-eQTLs were shared across three or more tissues at m > 0.9, 25.3% of FDR-matched cis-eQTLs were shared. The extensive tissue-specificity of trans-eQTLs was also supported by a hierarchical approach for FDR control[17], where we found no trans-eQTLs shared across more than one tissue (Extended Data Table 3). Our estimate of increased tissue specificity for trans-eQTLs agreed with the minimal sharing of trans effects reported in previous eQTL studies with fewer tissues[4,19], and greatly exceeds what would be expected on the basis of replication between tissues for cis-eQTLs of matched minor allele frequency (MAF) and effect size (Wilcoxon rank sum test; P ≤ 2.2×10−16 for all choices of replication FDR; Extended Data Fig. 8h). Given the greater tissue-specificity of trans-eQTLs, we note that heterogeneity in cellular composition of bulk tissue samples is one important confounder that may reduce power to detect trans-eQTLs, or even lead to false positive associations[6]. Despite the high tissue-specificity, we did observe a small number of tissue-shared trans-eQTLs, including rs7683255, which was moderately associated in trans with NUDT13 across most tested GTEx tissues with a consistent direction of effect (Extended Data Fig. 9a). We also found examples of trans-eQTLs shared across a subset of related tissues, such as an association between rs60413914 and RMDN3, a gene with increased expression levels in brain regions as compared to other tissue types, and for which the trans-eQTL had moderate effects in all tested brain regions but no strong effect in other tissues (Extended Data Fig. 9b, c).

Extended Data Table 3

trans-eVariant and trans-eGene discoveries with hierarchical FDr control

Tissue	No. of samples	No. of trans-eGenes	No. of trans-eVariants
Whole Blood	338	1	1
Skin – Sun Exposed (Lower leg)	302	2	3
Lung	278	2	2
Thyroid	278	2	2
Esophagus – Mucosa	241	3	3
Artery – Aorta	197	1	1
Skin – Not Sun Exposed (Suprapubic)	196	1	1
Heart – Left Ventricle	190	1	1
Testis	157	4	5
Colon – Sigmoid	124	1	1
Brain – Cortex	96	1	1
Brain – Putamen (Basal ganglia)	82	1	1
Total unique		20	22

Only tissues with non-zero discoveries are shown. The three-level hierarchical procedure (see Methods) performs FDR control across tissues. More specifically, it controls the FDR of eVariants, the average proportion of false variant-gene associations across all eVariants, and a weighted average of false tissue discoveries for the selected variant-gene pairs (weighted by the size of the eVariant and eGene sets). The procedure was applied after linkage disequilibrium pruning.

Extended Data Figure 9

Examples of trans-eQTLs shared across tissues

a, An example of a trans-eQTL (rs7683255–NUDT13) originally identified in sun-exposed skin (10%, P ≤ 1.1 ×10−10, indicated by asterisk) that has a global effect across tissues. The lines represent 95% confidence intervals of the effect sizes. A thicker line indicates that this variant–gene pair is called significant at P ≤ 0.05 in the corresponding tissue. b, An example of a trans-eQTL (rs60413914–RMDN3) that is genome-wide significant in putamen (basal ganglia) (10% FDR, P ≤1.2 ×10−13, indicated by asterisk) that has an effect in all five brain tissues tested but shows little effect in other tissues. c, Expression levels (RPKM) of RMDN3 in all donors across 44 tissues. Box plots depict the IQR, whiskers depict 1.5× IQR.

Multi-tissue cis-eQTL analyses have been shown to increase power by explicitly modelling sharing patterns across tissues[15,18,20]. We did not observe an improvement in power for trans-eQTL discovery, consistent with the limited sharing observed across tissues (Extended Data Table 3). However, we did observe improvements for cis-eQTL discovery, particularly among tissues with smaller sample sizes (Extended Data Fig. 10). To ensure that these findings did not depend on the modelling assumptions of Meta-Tissue, we analysed the P values for all genes and all tissues with treeQTL, which controls the FDR of eGene discoveries across tissues[17]. This procedure identified 17,411 cis-eGenes at 5% FDR, 2,314 fewer eGenes than with the single-tissue analysis. Although this analysis is more conservative overall than the tissue-by-tissue analysis, we observed an increase in the number of eGenes detected in the tissues with the smallest sample sizes, including several brain regions, as well as an increase in the average number of tissues in which an eGene was detected (from 7.8 for single-tissue analysis to 8.3; Extended Data Fig. 10). Additional cis-eGenes identified through meta-analysis were more likely to be significant as sample size increased compared to similar numbers of eGenes identified using a less stringent single- tissue FDR (Extended Data Fig. 10). This suggests that one strategy for increasing power in studies of inaccessible or sample-limited cell types would be to analyse them jointly with data from GTEx tissues.

Extended Data Figure 10

Sharing information across tissues for cis-eQTLs

a, The proportion of expressed genes for which cis-eGenes are discovered in single tissues (5% FDR; origin) and the multi-tissue meta-analysis (m > 0.9; arrow), stratified by the sample size of individual tissues. In the meta-analysis, cis-eQTL discoveries are made using Meta-Tissue to identify tissues where the posterior probability a given cis-eQTL effect exists (that is, the tissue’s m value for the variant) is >0.9. b, The proportion of expressed genes that had a cis-eQTL in the subsampled data (n = 70) is shown on the y-axis, and the actual sample size of the tissue is shown on the x-axis. The proportion is shown for the tissue-by-tissue approach (5% FDR; origin) and using Meta-Tissue (m >0.9; arrow). c, For each of the subsampled tissue data sets (n = 70), we identified the additional K discoveries that were made using Meta-Tissue but were not significant at the 5% FDR threshold in the tissue-by-tissue analysis; we then identified the K most significant cis-eQTLs in the tissue-by-tissue analysis with a q value greater than 5% representing the additional discoveries we would make by simply relaxing the FDR. We then compared these two sets of K cis-eQTLs to the list of significant cis-eQTLs found in the full tissue-by-tissue analysis by calculating the proportion of the K cis-eQTLs that were significant in the full analysis (y-axis). The tissues are ordered along the x-axis by increasing sample size in the actual data set. d, The proportion of annotated and expressed genes that were found to be eGenes using the tissue-by-tissue approach and the hierarchical selection procedure implemented by TreeQTL. e, The number of cis-eGene discoveries per tissue (y-axis) against sample size (x-axis). The number of discoveries for the tissue-by-tissue approach are represented by the origin of each segment, while the number of discoveries from the hierarchical selection procedure are shown as arrows. As with Meta-Tissue, the hierarchical procedure improves cis-eGene discovery for tissues with low sample sizes, albeit fewer tissues overall have the benefits of this effect. Furthermore, an outcome of this procedure is that for tissues with high sample sizes, reported numbers of cis-eGenes are more conservative than those observed in the tissue-by-tissue analyses or Meta-Tissue. f, Improvement of cis-eGene discovery by incorporating genomic annotations. For the 26 tissues for which we can relate cell-type specific chromHMM annotations, we identify cis-eGenes accounting for the variant-level genomic annotations and corresponding enrichment estimates using the Bayesian FDR control procedure described previously[80]. For each tissue, the number of cis-eGenes identified by the Bayesian procedure (arrow) is plotted against the tissue-by-tissue results (origin).

To characterize the biological properties of multi-tissue cis-eQTLs, we annotated discovered eVariants using chromatin state predictions from 128 cell types sampled by the Roadmap Epigenomics project[2]. eVariants were enriched in predicted promoter and enhancer states across all Roadmap cell types (Fig. 3a). However, the eVariants exhibited significantly greater enrichment in promoters and enhancers from matched tissues (Wilcoxon rank sum test, P ≤ 9.3 ×10−4, Extended Data Table 4), illustrating consistent patterns of tissue specificity for cis-regulatory elements and eQTLs (Fig. 3a). Furthermore, eQTL activity was significantly more likely to be shared across pairs of tissues when the eVariant overlaps the same chromatin state in both tissues (Wilcoxon rank sum test, P ≤ 5.0×10−5, Fig. 3b).

Figure 3

Functional characterization of cis-eQTLs

a, Enrichment (x-axis) of eVariants in cis-regulatory elements (CREs) across 128 Roadmap Epigenomics project cell types, for each GTEx discovery tissue (y-axis). Enrichment estimated by comparing to random MAF- and distance-matched variants. Stronger enrichment was observed in matched tissues (coloured dots) than in unmatched tissues (box plots). b, Proportion of eQTLs shared between two tissues (m > 0.9) if the eVariant overlaps the same Roadmap annotation in both tissues (y-axis) or different annotations (x-axis). Points represent the mean across all tissues, coloured by the discovery tissue. c, Enrichment of eVariants (y-axis) in tissue-matched enhancers (black) and promoters (grey) for the first four conditionally independent eQTLs discovered for each eGene (x-axis). d, Proportion of eVariants overlapping tissue-matched DNase I hypersensitive sites (DHS; y-axis) as a function of the probability that a variant is causal (x-axis), coloured by the eQTL discovery tissue. e, Normalized absolute eQTL effect size (x-axis) for each eVariant annotation class (y-axis). f, Median (line) and interquartile range (shading) of normalized absolute eQTL effect size (y-axis), as a function of the number of tissues in which the eGene is expressed (x-axis). Box plots depict the interquartile range (IQR), whiskers depict 1.5× IQR. OR, odds ratio.

Extended Data Table 4

GTEx tissue mapping with Epigenomics roadmap cell types

GTEx Tissue	Epigenomics Roadmap Cell Type
Adipose – Subcutaneous	Adipose Nuclei (E063)
Adipose – Visceral (Omentum)	Adipose Nuclei (E063)
Adrenal Gland	NA
Artery – Aorta	Aorta (E065)
Artery – Coronary	NA
Artery – Tibial	NA
Brain – Anterior cingulate cortex (BA24)	Brain Cingulate Gyrus (E069)
Brain – Caudate (basal ganglia)	Brain Anterior Caudate (E068)
Brain – Cerebellar Hemisphere	NA
Brain – Cerebellum	NA
Brain – Cortex	Brain Angular Gyrus (E067), Brain Inferior Temporal Lobe (E072), Brain Dorsolateral Prefrontal Cortex (E073)
Brain – Frontal Cortex (BA9)	Brain Inferior Temporal Lobe (E072), Brain – Dorsolateral Prefrontal Cortex (E073)
Brain – Hippocampus	Brain Hippocampus Middle (E071)
Brain – Hypothalamus	NA
Brain – Nucleus accumbens (basal ganglia)	NA
Brain – Putamen (basal ganglia)	NA
Breast – Mammary Tissue	Breast Myoepithelial Primary Cells (E027)
Cells – EBV-transformed lymphocytes	Lymphoblastoid Cells (E116)
Cells – Transformed fibroblasts	NA
Colon – Sigmoid	Sigmoid Colon (E106)
Colon – Transverse	Colonic Mucosa (E075), Colon Smooth Muscle (E076)
Esophagus – Gastroesophageal Junction	Esophagus (E079)
Esophagus – Mucosa	Esophagus (E079)
Esophagus – Muscularis	Esophagus (E079)
Heart – Atrial Appendage	Right Atrium (E104)
Heart – Left Ventricle	Left Ventricle (E095)
Liver	Liver (E066)
Lung	Lung (E096)
Muscle – Skeletal	Skeletal Muscle Male (E107), Skeletal Muscle Female (E108)
Nerve – Tibial	NA
Ovary	Ovary (E097)
Pancreas	Pancreas (E098)
Pituitary	NA
Prostate	NA
Skin – Not Sun Exposed (Suprapubic)	NA
Skin – Sun Exposed (Lower leg)	NA
Small Intestine – Terminal Ileum	Small Intestine (E109)
Spleen	Spleen (E113)
Stomach	Stomach Mucosa (E110), Stomach Smooth Muscle (E111)
Testis	NA
Thyroid	NA
Uterus	NA
Vagina	NA
Whole Blood	Primary mononuclear cells from peripheral blood (E062)

Integration of genomic annotations such as chromatin state has been demonstrated to improve power for eQTL discovery[8,21-23]. For 26 GTEx tissues matched with cell-type specific annotations from the Roadmap Epigenomics project, we applied a Bayesian hierarchical model for eQTL discovery by incorporating variant-level genomic annotations[24] that provided a substantial boost to discovery power. Inclusion of genomic annotations (enhancers, promoters and distance to the TSS) increased the total number of cis-eQTL discoveries by an average of 43% (or 1,200 genes) per tissue (Extended Data Fig. 10f), demonstrating the considerable advantage of integrating genomic annotations into eQTL mapping models.

Conditionally independent (secondary) cis-eQTL signals were located further from the TSS (median distance 50.1 kb from the TSS, compared to 28.9 kb for primary eQTLs; Wilcoxon rank sum test, P ≤ 2.2×10−16). However, similar to primary eVariant associations, secondary eVariants were enriched for chromosomal contact with target eGene promoters, as determined through Hi-C, compared to background variant–TSS pairs (Supplementary Information 6). This suggests that, despite their sequence-based distance from the TSS, primary and secondary eVariants are in close physical contact with their target gene promoters via chromatin looping interactions. While primary eVariants were significantly more enriched in promoters than enhancers, secondary associations exhibit increased enhancer enrichment, consistent with their increased distance from the TSS and tissue-restricted activity (Wilcoxon rank sum test, P ≤ 2.2×10−16, Fig. 3c).

To identify causal variants that are likely to underlie cis-eQTLs, we applied two computational fine-mapping strategies[25,26] (Supplementary Information 13 and Supplementary Figs 11, 18). First, we identified 90% credible sets (that is, the collection of variants with 90% probability of containing all causal variants) for each eGene in each tissue using CAVIAR[25]. Across all tissues, the mean credible set size was 29 variants (per tissue means ranged from 25 to 31). Second, we estimated the probability that each eVariant is a causal variant using CaVEMaN[26]. Across tissues, between 3.5% and 11.7% of top eVariants were predicted to be causal variants (causal probability P > 0.8). Consistent with variants with high causal probabilities being functional regulatory variants (as opposed to linkage disequilibrium proxies), 24.3% of eVariants with causal probabilities in the top tenth percentile (0.77