Introduction

Individual differences in RNA expression may result from variation in environmental exposures, stochastic variation, age, sex and genotype differences [1] and thereby may also be involved in the widely observed contribution of DNA, age and sex to the etiology of complex diseases [2,3,4,5]. Quantifying human transcriptomic heritability (h2) is of interest for transcriptomic and genomic studies. For example, one possible reason for the association of gene expression with a certain phenotype is that DNA variants influence the amount of gene expression through expression quantitative trait loci (eQTLs). Gene expression can have substantial h2 not explained by eQTLs [6], and therefore the h2 of gene expression and corresponding eQTL findings can be considered complementary to transcriptomic and genomic studies.

The h2 of whole blood gene expression has been established previously for genome-wide transcriptomic data generated by microarray technology and by RNA-Seq. For array technology, Wright et al. [7] found a mean gene expression h2 of 0.10 (SD = 0.14, N = 2752, 18,392 genes) from modeling data assessed in mono- and dizygotic twin pairs. Local identity-by-descent (IBD) analyses, which provide an estimate of variance of gene expression explained by genetic relatedness, resulted in a mean local h2 of 0.03, explaining 23% of the total heritability. A population-based study (N = 2765) by Lloyd-Jones et al. [6] found mean estimates of h2 of microarray-based gene expression of 0.19 with a local h2 of 0.06, resulting in a mean proportion of genetic variance explained by all eQTLs of 31%. Viñuela et al. [8] estimated a mean of blood-derived h2 of 0.23 in a twin-based sample (N = 855) for RNA-Seq-based data.

Here, we analyzed RNA-Seq (Illumina Hiseq2000) data from 459 MZ twin pairs, 150 DZ twin pairs, 24 relatives of twin pairs, and 255 unrelated participants, leading to a total dataset with genotype and expression data of 1497 adult participants (998 females) from the Netherlands Twin Register [9,10,11] (Supplementary Table 1). Our RNA-Seq dataset contained a large (>60%) sample overlap with the microarray-based dataset from Wright et al. [7], allowing a reliable comparison between microarray and RNA-Seq-based h2 estimates. We first estimated twin-based h2 of gene expression by making use of the classical twin design [12]. Next, we simultaneously estimated the variance explained by a genetic relationship matrix (GRM)-containing SNPs in a 250 kilobase (kb) cis-window of the gene (h2cis), and the variance explained by a second GRM including only closely related individuals (h2res) using GCTA [13]. Together, h2cis and h2res constitute the total heritability (h2total). We performed cis- and trans-eQTL analyses of the same samples and compared h2 and eQTL findings to test for consistency.

Methods

Participants

RNA samples were obtained from 1497 participants from the Netherlands Twin Registry (NTR) and included 459 complete MZ twin pairs and 150 complete DZ twin pairs, 24 relatives of twin pairs, and 255 unrelated participants. NTR is a longitudinal cohort study of twins and their families [10, 14]. The age of the participants ranged from 17.6 to 79.6 years old (mean = 36.7, SD = 14.0), 67% of the sample was female. The data used for this study largely overlap (60%) with those used in an earlier study [7]. See Supplementary Table 1 for a description of the samples.

RNA extraction and sequencing

Venous samples were drawn in the morning after an overnight fast. Heparinized whole blood samples were transferred within 20 min of sampling into PAXgene Blood RNA tubes (Qiagen, Valencia, CA, USA) and stored at −20 °C. Total RNA from whole blood was depleted of globin transcripts using the Ambion GLOBINclear kit and subsequently processed for sequencing using the Illumina TruSeq version 2 library preparation kit. Paired-end sequencing of 2 × 50-bp reads was performed using the Illumina HiSeq 2000 platform, pooling ten samples per lane and aiming for >15 million read pairs per sample. Adapters were identified and clipped, and low-quality read ends were trimmed (min length 25, min quality 20 Read alignment was performed using STAR 2.3.0e against the Genome of the Netherlands (GoNL) reference panel [15]. Expression was using Ensembl v.71 annotation (which corresponds to GENCODE v.16). Overlapping exons (on either of the two strands) were merged into meta-exons, and expression was quantified for the whole meta-exon, resulting in base counts per exon or meta-exon. Gene expression, as base count per gene, was calculated as the sum of the expression values for all exons of each gene (excluding meta-exons). This pipeline is explained in detail in Zhernakova et al. [16].

Gene expression values per gene were ranked and mapped to a normal distribution with mean 0 and SD 1, after which values were corrected for sex, age, and cell counts: monocyte, lymphocyte, eosinophil, basophil, neutrophil and red blood cell counts, and 27 measurement batches. We then performed a principal-component analysis, and values were corrected on only the first principal component, which was not heritable based on comparison of goodness of fit between different twin-based structural equation models (p = 0.74). This PC explained 16.5% of the variance in the data (see Supplementary Figs. 1 and 2). Finally, expression values were centered and subsequently scaled by dividing these values by their respective standard deviations. Analyses of both classical twin modeling-based and identity-by-state (IBS)-based h2 were based on this final dataset.

Genotype data

Within the NTR, genotype information is available for 15,111 individuals for four different genotyping arrays (Affymetrix 6.0 (N = 11,781), Affymetrix Perlegen 5.0 (N = 1265), Illumina 660 (N = 1439) and Illumina Omni Express 1 M (N = 257), as well as sequence data from the Netherlands reference genome project GONL (BGI full sequence at 12 × (N = 368; [17]). Samples were removed if they had a genotype call rate below 90%, heterozygosity fell outside the range of −0.075 to 0.075, gender and IBS status mismatch occurred, or if the Mendelian error rate was larger than 5 standard deviations from the mean of all samples and for samples measured on Affymetrix 6.0 when the contrast quality control value was smaller than 0.40. Quality control of the SNPs was done for each platform separately, with SNPs being removed when they could only be aligned to the forward strand of the reference panel, the allele frequencies differed more than 10% with the reference allele, minor allele frequency (MAF) was below 0.005, Hardy–Weinberg equilibrium (HWE) test p < 10−12, and a genotype call rate of <0.95. The data of the different genotyping methods, except GONL sequence individuals, were subsequently merged into a single dataset. The missing SNP genotypes between each platform were imputed to the GONL reference data. Filtering of the imputed dataset included the removal of SNPs which were significantly associated with a single genotyping platform, if the allele frequency difference differed more than 10% with the GONL reference set, HWE p < 10−5, Mendelian error rate >mean + 5 SD and if the imputation quality (R2) was below 0.90. After filtering the GONL samples were readded to a cross-platform imputed dataset that includes 1,261,818 SNPs. We did not perform additional genotype quality control after subsetting the individuals with RNAseq.

Statistical methods

We employed two methods for estimating total heritability, a classical twin modeling-based approach, and an IBS-based approach. The classical twin modeling approach requires information on relatedness between subjects, while the IBS-based approach quantifies the genetic relatedness between subjects based on their genome-wide genetic data. For both methods, only autosomal genes were considered.

Classical twin modeling

Based on the resemblance of MZ and DZ twin pairs, the variance of the expression of each gene can be decomposed into additive genetic, common (or shared) environmental and unique environmental variance. Classical twin-based modeling was done in the structural equation modeling (SEM) R-package OpenMx [18, 19]. Models were fitted to decompose the variance of gene expression due to additive genetic (A), shared environmental (C), and unique environmental (E) effects. Additive genetic effects combine all the effects of genetic variants influencing gene expression. Shared environmental variance represents the proportion of variance explained by effects are shared by both members of a twin pair. Unique environmental variance results from environmental effects that are not shared by twins. We used a standard ACE model assuming dizygotic twins have an average IBD sharing of 0.5 across the genome, and monozygotic twins share an IBD of 1. Parameters were estimated by maximum likelihood (ML). We restricted the estimates to be positive. The results of these analyses provide the ML estimates of variance components. The comparison of ACE and CE models gave an estimate for the significance of the A-component.

IBS-based analysis

Techniques to quantify genetic similarity of ‘unrelated’ individuals who are genotyped for a large number of SNPs across the entire genome have been developed to estimate h2 due to SNPs. In software packages such as GCTA [13] the relatedness among individuals based on measured SNPs can be combined with known genetic relatedness in relatives in a two-variance component linear mixed model [20], in which simultaneous estimation of SNP heritability and total h2 is feasible. We used NTR genotype information for the 1497 individuals for whom both expression- and genotype data were available. A GRM was created for each cis-window of a gene, defined by the coordinates of a gene with 250 kb flanking area on each side. This GRM is referred to as the cis-GRM and represents the variance explained by all measured SNPs (and SNPs tagged by these measured SNPs) in a cis-window around the gene of 250 kb. Variance explained by this cis-GRM is referred to as cis-h2 (h2cis). The mean number of SNPs in cis-GRMs was 223.2 (SD = 119.9) (see Supplementary Fig. 3). A second GRM including closely-related individuals (that is, a genetic correlation > 0.05) was created for the autosomes. This GRM is referred to as the residual GRM. This GRM had all off-diagonal elements below 0.05 set to 0, to remove distant relatedness from the matrix. Variance explained by this residual GRM is referred to as residual h2 (h2res). The sum of h2cis and h2res is h2total (h2total = h2cis + h2res). As the cis-GRM is based on a limited number of SNPs there is substantial power to detect the genetic effects in cis [21], the second GRM will absorb all genetic variance not explained by the SNPs in the cis-window, or in high LD with SNPs in the cis-window. Note that due to the presence of a large number of related individuals this GWAS will capture genetic variance tagged by substantial IBD sharing and thus the sum of the two effects will be approximately equal to the heritability estimated in a twin study.

$${\mathrm{cov}}\left( {\mathrm{expression}} \right)_{n \ast n} = {\mathrm{GRM}}_{n \ast n}^{\mathrm{IBS}}{\mathrm{\Theta }}\sigma _{cis - {\mathrm{SNPs}}}^2 + {\mathrm{GRM}}_{n \ast n}^{\mathrm{IBS} > 0.05}{\mathrm{\Theta }}\sigma _{\mathrm{SNPs}}^2 + {I}_{n \ast n}{\mathrm{\Theta }}\sigma _{\mathrm{e}}^2$$

A total of 52,844 genes were analyzed, and subsequently filtered for being protein coding, having read counts above zero in at least 85% samples in each zygosity group (e.g. expressed in in at least 780 MZ twins and 255 DZ twins), a median expression count above 10, and more than 20 SNPs in the cis-window, resulting in an analysis of 11,353 genes (Supplementary Table 2).

Annotation and enrichment

We employed multiple annotation steps to interpret h2 estimates. We tested whether h2 was correlated with gene expression level, gene length, GC content or several loss-of-function scores obtained from Lek et al. [22] using linear regression. Gene locations and lengths were downloaded on 2017-12-28 using the Biomart community portal [23] using data from Ensembl [24]. We investigated genes for which the expression correlates highly (>0.8) with other genes.

In addition, we tested whether heritable genes are overrepresented in the canonical gene pathways from the molecular signature database (MSIGDB):: KEGG, REACTOME, BIOCARTA pathways downloaded from http://software.broadinstitute.org/gsea/downloads.jsp (c2.cp.v6.1) and genes identified in GWAS for immune diseases, mental or behavioral disorders, cardiovascular diseases, or cancer (extracted from the GWAS catalog [25] (as of September 2018) with the search terms: ‘immune system disease’, ‘mental or behavioral disorder’, ‘cardiovascular disease’ or ‘cancer’, respectively). For each pathway/gene group, a Wilcoxon test was performed between the median heritability of the genes in the pathway and the median heritability of the genes outside the pathway. Comparisons were made with previous analyses of h2 of gene expression from a recent whole blood-based RNA-Seq study by Battle et al. [26] and the GTEx project [27], which published heritability estimates for gene expression in adipose tissue (subcutaneous), tibial artery, heart (left ventricle), lung, muscle (skeletal), tibial nerve, skin (sun-exposed), thyroid, and whole blood.

eQTL analysis

The 1497 gene-level RNA samples with the same preprocessing as described above which were also used for heritability analysis, were also used for cis- and trans-eQTL analysis. For this analysis, the same cross-platform imputed dataset that includes 1,261,818 SNPs as described above, was filtered at MAF >0.01 and HWE <1 × 10−3, resulting in 1,239,670 SNPs.

For cis-eQTL analysis all associations between DNA variants and genes at distance <250 kb were computed, for trans-eQTL analysis all SNP - gene pairs at distance >250 kb. eQTL effects were detected with a mixed linear model approach using fastGWA a implemented in GCTA (https://cnsgenomics.com/software/gcta/#fastGWA, https://www.biorxiv.org/content/10.1101/598110v1).

For fastGWA, first a GRM is built using the –make-grm option in GCTA. Then, a sparse GRM is built using the option –make-bK-sparse 0.05. fastGWA is then run using this sparse GRM, with expression level as the dependent variable and SNP genotype values as independent variable. Correction for multiple testing was done using FDR, for cis- and trans-eQTL analysis separately, resulting in a P-value threshold of 1 × 10−5 for cis-eQTLs, and 1.5 × 10−7 for trans-eQTLs. We are aware that FDR used like this may result in more false positives than 5%, however, we are merely interested in the overlap between heritability and eQTL analysis and do not draw any conclusions on the amount of identified eQTLs.

Simulation

To gain insight in the performance of our models, we simulated twin-based phenotype data with prespecified heritabilities and tested whether our estimations were in accordance. We ran OpenMX and GCTA models using simulated phenotypes and real genotype data, with different values for variance components (see Supplementary Figs. 6 and 7).

Results

The gene-level whole blood RNA-Seq data from 1497 participants, including 459 MZ twin pairs and 150 DZ twin pairs, originally contained 52,844 genes (Supplementary Table 2). After filtering (see “Methods”), 11,409 protein coding genes that were expressed in at least 85% of the samples were left for analysis. In these genes, twin-based h2 for each gene was estimated based on the classical twin design (see Methods). Next, IBS-based methods were applied to compute SNP heritability.

For the twin-based analyses, a genetic structural equation model [28] that included additive genetic, common, and unique environmental factors was fitted to the expression data from each gene. The mean h2 (the standardized estimate of the contribution of additive genetic factors) was estimated at .20, the standardized mean contribution of shared environment at 0.05, and the standardized mean unshared environment at 0.74. There was a considerable spread in estimates, with estimates ranging between 0 and 1, except for the common environment component which had a maximum of 0.51 (Fig. 1 and Table 1).

Fig. 1
figure 1

Histograms of estimates of heritability of gene expression. a Classical twin modeling-based estimates of each gene. Based on the resemblance of MZ and DZ twin pairs, the variance of the expression of each gene was decomposed into additive genetic, common (or shared) environmental, and unique environmental variance. A standard ACE model was used, assuming dizygotic twins have an average IBD sharing of 0.5 across the genome, and monozygotic twins share an IBD of 1. Parameters were estimated by maximum likelihood (ML). b IBS-based h2cis estimates of each gene, i.e. variance explained by a GRM created for each cis-window of a gene, defined by the coordinates of a gene with 250 kilobase (kb) flanking area on each side. c IBS-based h2res estimates of each gene, i.e. variance explained by a GRM including closely related individuals for the autosomes with all off-diagonal elements below 0.05 set to 0

Full size image
Fig. 2
figure 2

Spearman correlations between h2cis and estimates of GTEx per-tissue heritabilities [27]. a Tissue-specific heritabilities. b Tissue-wide heritabilities

Full size image
Table 1 Estimates for classical twin modeling-based and IBS-based heritability of gene expression
Full size table

Since estimates for the contribution of common environment on average were low, we proceeded with an IBS-based approach that did not take into account common environment. Gene expression is controlled by both local (cis) eQTLs and non-local (trans) eQTLs. Therefore, when studying the genetic component of expression in terms of heritability, we also make a distinction between the genetic component that is close to the gene (h2cis), and the genetic component that is not close to the gene (h2res). h2cis consisted of the variance explained by SNPs in a 250-kb cis-window of the gene and h2res the variance explained by genome-wide close relatedness. From the IBS approach we obtained heritability estimates for 11,353 out of 11,409 genes. The correlation between IBS-based total h2 (h2cis + h2res) and classical twin modeling-based h2 was 0.98 (p < 10−308, Supplementary Fig. 4). The IBS-based approach resulted in a mean h2res of 0.20 and a mean h2cis of .06 (Fig. 1 and Table 1), summing up to a h2total of 0.26. We found 721 genes to have a Bonferroni-corrected significant h2cis (p < 4.40 × 10−6), all of which had a h2cis larger than 0.16. We found 5636 genes to have a Bonferroni-corrected significant h2res (p < 4.40 × 10−6), all of which had a h2res larger than 0.01 (Supplementary Table 3, Supplementary Fig. 5). The mean contribution of h2cis to the total h2 (mean h2cis/(mean h2cis + mean h2res)) was 20%, with a range from 0 to 100%. The correlation between h2res and h2cis was 0.06 (p = 6.36 × 10−10). We found good performance of our models, slightly underestimating local heritability (see Supplementary Figs. 6 and 7).

Correlation of gene expression h 2 between tissues

A whole blood-based RNA-Seq study by Battle et al. [26] published heritability estimates due to regulatory genetic variation (N = 922). Estimates of heritability due to local genetic variation from this study correlated 0.81 with h2cis from our study (p < 10−308, Spearman correlation, Supplementary Fig. 8). eQTL studies have shown that gene expression in different tissues is regulated by DNA by partially overlapping, and partially unique QTLs (https://science.sciencemag.org/content/348/6235/648). To study the extent to which the total local genetic component of gene expression is similar between tissues, we looked at h2 estimates of gene expression in The Genotype-Tissue Expression (GTEx). The GTEx consortium, (N = 422) [29] reported h2 estimates for gene expression from RNA-Seq in adipose tissue (subcutaneous), tibial artery, heart (left ventricle), lung, muscle (skeletal), tibial nerve, skin (sun-exposed), thyroid, and whole blood, estimating local h2 both unique for a tissue (tissue specific) and heritability shared between tissues (tissue wide). Our estimates of h2cis showed significant correlations with local h2 in every measured tissue, in both tissue-wide (ρ > 0.12) and tissue-specific heritability (ρ > 0.03) estimates. The strongest correlation was found between h2cis and tissue-wide heritability of whole blood-derived gene expression (ρ = 0.24, p = 1.99 × 10−144, Fig. 2, Supplementary Table 4). Heritability estimates reported in the GTEx study correlated <.44 between tissues in tissue-wide estimates and <0.19 in tissue-specific estimates.

Fig. 3
figure 3

Heritability estimates for genes divided into groups based on gene expression levels per decile. a Total heritability estimates from Wright et al. [7] b IBS-based h2total estimates. c IBS-based h2cis estimates. d IBS-based h2res estimates

Full size image

Correcting h2cis for median read count only showed marginal effects on correlations. There were 614 genes that were highly expressed in all tissues in GTEx (above the 8th decile). These genes showed decreased h2cis (p = 5.58 × 10−12, mean = 0.04, median = 0.02) and increased h2res (p = 7.77 × 10−14, mean = 0.23, median = 0.22) compared to the full set of tested genes.

Pathway analyses

To annotate gene expression heritability, we studied if certain gene pathways have higher average heritability than expected. In order to so we tested for enrichment of h2res and h2cis in canonical gene pathways covering a broad range of biological pathways that are well curated (KEGG, REACTOME, BIOCARTA). The expression of genes identified in GWAS is likely to be under genetic control: to test if this is the case for gene expression in blood, we also tested if h2res and h2cis is enriched in genes identified in genome-wide association studies (GWAS) for immune diseases, mental and behavioral disorders, cardiovascular diseases, or cancer (extracted from the GWAS catalog [25]) to cover GWAS findings for a broad range of diseases. Enrichment analysis were performed before and after correcting h2 for median gene expression per gene.

We observed significantly higher h2res (false discovery rate (FDR) <0.05) in 343 canonical pathways (top hit: KEGG cytokine-cytokine interaction, p = 4.22 × 10−15), and in genes identified in GWAS for immune diseases (p = 4.99 × 10−13), mental disorders (p = 1.79 × 10−7), cancer (p = 5.53 × 10−5) and cardiovascular diseases (p = 1.06 × 10−9). After correction for mean gene expression, h2res was significantly higher in 125 canonical pathways (top hit: KEGG cytokine-cytokine interaction, p = 6.66 × 10−15) and in genes identified in GWAS for immune diseases (p = 2.33 × 10−10), mental disorders (p = 3.63 × 10−10), cancer (p = 7.10 × 10−6) and cardiovascular diseases (p = 9.69 × 10−10).

We found significantly higher h2cis (FDR <0.05) in 6 canonical pathways (top hit: KEGG lysosome, p = 1.41 × 10−7), and in genes identified in GWAS for immune diseases (p = 6.71 × 10−4), mental disorders (p = 7.13 × 10−4), cancer (p = 1.09 × 10−7) and cardiovascular diseases (p = 3.24 × 10−5). After correction for mean gene expression, h2res was significantly higher in 10 canonical pathways (top hit: KEGG lysosome, p = 5.79 × 10−9) and in genes identified in GWAS for immune diseases (p = 9.13 × 10−5), for mental disorders (p = 2.36 × 10−3), cancer (p = 3.15 × 10−7) and cardiovascular diseases (p = 1.16 × 10−5) (see Table 2 and Supplementary Table 5).

Table 2 Gene pathway enrichment of heritable genes in certain gene pathways (KEGG, REACTOME, BIOCARTA) and genes identified in GWAS for immune diseases, mental or behavioral disorders, cardiovascular diseases, or cancer (extracted from the GWAS catalog [25])
Full size table

Gene expression h 2 correlations

In order to identify physiological gene properties that are correlated with gene expression heritability, we evaluated the correlation between heritability and gene expression level, gene length, GC content and several loss-of-function scores, see Table 3). We found a significant association of h2res and h2cis with median read count (P = 6.81 × 10−276, P = 2.15 × 10−2, respectively) and GC content of a gene (P < 1.80 × 10−115, P = 4.02 × 10−26, respectively). After correcting for median read count, GC content was still significantly correlated with h2res and h2cis (P = 2.27 × 10−3, P = 6.11 × 10−28, respectively). The length of a gene was significantly correlated with h2res, with longer genes having a slightly higher h2res (P = 4.94 × 10-8). Gene length did not influence h2cis (P = .50). A high intolerance to LoF or high probability of loss-of-function (pLI) did not significantly influence h2 estimates.

Table 3 Predictors of h2
Full size table

Relation of h 2 cis with strength of cis-eQTLs

Gene expression can have substantial h2 not explained by eQTLs [6]. In order to study the overlap between heritability and eQTL results, we performed eQTL analysis in the same sample (see Methods) and found 5249 genes with a significant cis-eQTL (p-value threshold 1.5 × 10−7 for a FDR of 5%). In addition, we found a significant association between h2cis and the beta of the corresponding top cis-eQTLs (ρ = 0.7644, p < 10−308). We also tested the correlation between h2cis and the presence of cis-eQTLs in results from Zhernakova et al. [16], who performed RNA-Seq-based eQTL analysis in an independent sample (N = 2116 unrelated adults). There was a strong correlation between h2cis and the Z-score of the strongest eQTL (ρ = 0.75, p < 10−308) (Supplementary Fig. 9).

Relation of h 2 res with strength of trans-eQTLs

We performed eQTL analysis in the same sample for SNPs outside a cis-window of 250 kb around each gene, resulted in 2433 genes with a significant trans-eQTL (p-value threshold 1.5 × 10−7 for a FDR of 5%). Estimates of h2res correlated with the beta of the corresponding top trans-eQTLs (ρ = .04, p = 1.89 × 10−3). The mean p-value of trans-eQTL of genes with h2res > 0.4 was 5.51 × 10−9 (median = 1.44 × 10−9).

Comparing RNA-Seq and micro-arrays: relation of mean expression levels with heritability

The difference in heritability estimates between RNAseq and array data results, especially with gene expression intensity as a factor, is a valuable comparison to distinguish noise or bias that is inherent to a specific technique. As RNAseq is a more expensive method, it is informative to know if and in which situations this technique offers advantages to answer a particular research question. From our original results on 52,844 genes, we selected the 12,070 genes that were measured both with Affymetrix expression arrays in the study by Wright et al. [7] and with RNA-Seq with read counts above zero in at least 85% of samples in each zygosity group. In contrast to previous analyses, we did not apply any additional filtering (i.e. minimum median read count per gene, protein coding genes only, minimum of SNPs per GRM to allow for an unbiased comparison. In the RNA-Seq data, the Spearman correlation of mean expression with h2total was −0.13 after correcting for covariates (see Methods) (p = 7.49 × 10−45). In the array-based results from Wright et al. [7] the Spearman correlation between mean expression and h2 was 0.28 (p = 2.08 × 10−219) after correcting for covariates.

Using this selection of 12,070 genes, mean h2res, h2cis, h2total, and h2total results from Wright et al. [7] were 0.19, 0.06, 0.25, and 0.14, respectively. When we partitioned mean expression levels into 10 deciles for both microarray and RNA-Seq data, we saw that h2res estimates were also higher in RNA-Seq data for almost all deciles (see Fig. 3 and Supplementary Table 6), with the difference being largest in the lowest deciles (p = 1.06 × 10−72). In the highest decile of expression, the h2 of array-based expression was higher (p = 2.00 × 10−6). This suggests that the resolution of RNA-Seq is better able to capture variation in low to moderately expressed genes. Genes measured by both array and RNA-Seq were mostly in the same or nearest decile (see Supplementary Fig. 10). Estimates of h2cis showed a slight negative correlation with median expression level in the RNA-Seq data eQTLs (ρ = −0.04, p < 5.79 × 10−5).

Discussion

The present study estimated the h2 of gene expression in RNA-Seq-based expression data by making use of the different genetic relatedness of mono- and dizygotic twins and an IBS approach. The mean of total gene expression h2 (0.26) was substantially higher than found with the microarray-based study by Wright et al. with largely overlapping RNA samples. This was also the case with a direct comparison of 12,070 genes that were measured both with Affymetrix expression arrays in the study by Wright et al. [7] and with RNA-Seq, where mean RNA-Seq-based h2total was 0.25 and mean microarray-based h2total was 0.14 (p < 10−308). Heritability estimates in RNA-Seq did not increase with gene expression level, as opposed to the results from microarray data. This suggests RNA-Seq measurements are less noisy, in particular in genes with low expression, as compared to microarrays measures.

We estimated h2 of gene expression at 0.20 for mean h2res, and at 0.06 for mean h2cis. This resulted in a relative contribution of h2cis to h2total of 20%. This is in line with earlier findings by Wright et al. [7] (relative contribution of local IBS-driven h2 of 23%) and Lloyd-Jones et al. [6] (proportion of h2 explained by cis-eQTL of 0.31). Since local variants in the cis-window of a gene only explain 20% on average of the total h2, and the number of genes for which the majority of heritability stems from genetic relatedness in the cis-window is very low, our findings indicate that loci outside the cis-window of a gene or rare local variants explain a significant proportion of total h2 of gene expression. This is strengthened by our finding that the strength of a trans-eQTL is correlated with h2res.

The h2cis estimates correlated strongly with estimates from an independent sample (ρ = 0.81, whole blood, N = 922) [26], but much less so in the smaller GTEx dataset [27] (ρ = 0.24, whole blood, N = 449). This shows that reasonably large sample sizes are needed for accurate h2 estimation. However, even with the small sample size used in GTEx we found significant correlations between h2cis in whole blood and h2cis in the other tissues, which is in line with the finding that cis-eQTLs are partially shared between tissues [29].

We found a significantly higher h2res in 125 canonical pathways, with the strongest enrichment for genes in the KEGG cytokine-cytokine interaction pathway, and many other immune system pathways (including KEGG Innate Immune system (p = 8 × 10−12), REACTOME interferon signaling (p = 2 × 10−8) and KEGG natural killer cell cytotoxicity (p = 2 × 10−7)) [30].

Both h2cis and h2res were higher in all genes identified in GWAS for immune diseases, mental disorders, and cardiovascular diseases, although interestingly the enrichment was much stronger for h2res. This indicates that for genes associated with a disease through GWAS, expression in blood is not only locally regulated but also enriched with genome-wide SNP signal as reflected in the high average h2res, suggesting that SNPs found in GWAS are influencing expression of genes outside the cis-window of the gene. A high intolerance to loss-of-function (LoF) or high probability of loss-of-function (pLI) did not significantly influence h2 estimates, suggesting that increased mutational load in a gene increases the genetic variation as much as it influences variation in gene expression and therefore does not influence h2.

The comparison of mean h2 of genes across different RNA measurement techniques can be viewed as a proxy for the comparison of the measurement error between techniques. An increase in measurement error always implies a decrease in h2 (because measurement error introduces random divergence within twin pairs). Consider for example the h2 of probes conditioned on the median expression level. When comparing the h2res estimates obtained based on RNA-Seq with those obtained based on Affymetrix expression arrays [7], it becomes apparent that micro-array-based estimates of gene expression heritability are associated with gene expression levels. The fact that differential measurement error conditional on expression level plays less of a role for RNA-Seq data ensures that variation across genes reflects biological signal. The slight negative correlation we found of h2cis with median expression level in the RNA-Seq data eQTLs (ρ = −0.04, p < 5.79 × 10−5) is counterintuitive, difficult to interpret and presumably not meaningful.

If, or when, researchers eventually examine the h2 of RNA expression levels in single cells, or nuclei, an inspection of the relationship between median expression levels and h2 can be used to detect this source of differential measurement error.

In summary, this study shows possible advantages of h2res-informed trans-eQTL discovery, reproducibility of h2cis, and the benefits of using RNA-Seq for estimating heritability of low-expressed genes.