A mega-analysis of expression quantitative trait loci (eQTL) provides insight into the regulatory architecture of gene expression variation in liver

Genome-wide association studies (GWAS) have identified numerous genetic variants in the human genome associated with diseases and traits. Nevertheless, for most loci the causative variant is still unknown. Expression quantitative trait loci (eQTL) in disease relevant tissues is an excellent approach to correlate genetic association with gene expression. While liver is the primary site of gene transcription for two pathways relevant to age-related macular degeneration (AMD), namely the complement system and cholesterol metabolism, we explored the contribution of AMD associated variants to modulate liver gene expression. We extracted publicly available data and computed the largest eQTL data set for liver tissue to date. Genotypes and expression data from all studies underwent rigorous quality control. Subsequently, Matrix eQTL was used to identify significant local eQTL. In total, liver samples from 588 individuals revealed 202,489 significant eQTL variants affecting 1,959 genes (Q-Value < 0.001). In addition, a further 101 independent eQTL signals were identified in 93 of the 1,959 eQTL genes. Importantly, our results independently reinforce the notion that high density lipoprotein metabolism plays a role in AMD pathogenesis. Taken together, our study generated a first comprehensive map reflecting the genetic regulatory landscape of gene expression in liver.

regions are often intergenic but nevertheless can have an effect such as recruiting transcription factors, which in turn can influence expression of nearby genes 11 . In general, such loci potentially harbour regulatory sequences in cis or trans to the gene regulated by the associated genetic variant.
Correlating the allele count at a variant locus and the expression of nearby genes in a given tissue can bridge the gap between the observed genetic association and understanding the mechanisms responsible for disease risk by defining an expression quantitative trait locus (eQTL) 12 . In recent years, thousands of eQTL were identified in multiple tissues by genome-and transcriptome-wide approaches 13 . Disease-associated genetic markers that represent a significant eQTL for a nearby gene can thus easily be identified. For AMD, so far only a single eQTL (rs79037040) affecting the expression of the tumor necrosis factor receptor superfamily, member 10a (TNFRSF10A) in white blood cells was reported to be associated with disease risk 14 . The lack of additional eQTL involved in AMD pathology can possibly be attributed to the observation that many eQTL studies are greatly underpowered 15,16 . In addition, although around 50% of known eQTL are common to several tissues 13 , many eQTL are likely to be specific for a given tissue or cell type.
The primary site of disease in AMD is the retinal tissue complex consisting of the retinal pigment epithelium (RPE), Bruch's membrane and the choriocapillaris. The function of the liver is fundamentally different from the retina; thus the liver likely will react differently to environmental influences than retinal tissue. Furthermore, eQTL in liver might behave differently in retinal cells. However, it is challenging to sample a large number of human retinae and, as a consequence, no eQTL data from one of these cell types have been reported to date. Thus, we aimed at performing eQTL analysis in a surrogate tissue which expresses several genes of interest in loci associated with AMD, with the assumption that a polymorphism could have similar effects on gene expression in the surrogate tissue as in the retina. We selected liver as surrogate tissue since it is the main tissue for expression of genes of the complement system and of HDL metabolism. Moreover, gene products (e.g. proteins) of complement and of HDL metabolism expressed by the liver are frequently secreted into circulation where they exert various biological activities, and which could consequently influence AMD through its systemic effect in the choriocapillaris. With this rational we anticipated that investigating eQTL of these genes in liver could reveal important mechanistic insights into the association of these loci with AMD.
Several previous studies have published eQTL from liver tissue using different genotyping and expression profiling platforms [17][18][19][20] . Raw or curated data files of these studies are publicly available. In the present study, we have jointly analysed the data from the four independent liver eQTL resources by state-of-the-art methods, subsequent to rigorous quality control. In addition, the results were compared to published GWAS data for AMD risk variants. We show that a common, AMD associated deletion of the complement factor H related 1 and 3 genes (CFHR1/3) results in a markedly reduced expression of both genes in the liver. Furthermore, we show that two AMD risk variants are significant eQTL in liver affecting the expression of two genes involved in HDL metabolism.

Results
Data preparation. The main objective of this study was to identify significant cis-eQTL in liver tissue as part of our long-term goal to understand the functional consequences of genetic variants associated with complex diseases such as AMD. To this end, individual datasets publically available were merged although each one used distinct platforms to call genotypes and to measure gene expression (Table 1). Consequently, stringent quality control measures were applied to compile a data set of high quality genotypes and gene expression values comparable across studies. Altogether, the study comprised 6,256,941 imputed variants and expression values of 24,123 genes in 588 samples of European descent. eQTL Analysis. First, we performed eQTL calculations for each of the four studies individually 13,[17][18][19] . Local eQTL were calculated by including all variants on the same chromosome that are located within 1,000,000 base pairs (1 Mbp) up-or downstream of the transcription start site or polyadenylation site of a gene locus, respectively. Next, mixed effects models were used to perform a meta-analysis by including the effect sizes and standard errors obtained from each study separately. In order to account for multiple testing, we controlled the false discovery rate (FDR) to be smaller than 0.001 21 . At this threshold, 101,148 eQTL variants and 1,313 genes differentially regulated by the eQTL were identified ( Table 1).
As meta-analysing data can result in a loss of statistical power [22][23][24] , we additionally performed a mega-analysis by directly estimating eQTL in the entire dataset comprising all four studies. The mega-analysis yielded 202,489 statistically significant eQTL variants affecting the expression of 1,959 genes while controlling the FDR to be less than 0.001 (Fig. 1, Table 1 and Supplementary Table S1). Compared to the results from the meta-analysis, the mega-analysis provided a two fold increase in the number of eQTL variants and a 1.5 fold increase in the number of differentially regulated genes. Of note, however, both mega-and meta-analysis discovered more significant results than any of the four individual studies alone (Table 1). Only 38.5 to 60.9% of the significant single study eQTL genes could be replicated in the meta-analysis. The GTEx study had the lowest replication rate, possibly due to its relatively small sample size (N = 83). The overlap of single study results and the mega-analysis is on average 19% higher (53.5 to 80.15%) than the overlap observed in the meta-analysis. As the mega-analysis reproduced 95.96% of the meta-analysis eQTL and detected many signals beyond, we decided to rely on the data of the mega-analysis for further calculations although this may represent a slight overestimation of eQTL derived from the available data set.
We next aimed to identify independent eQTL variants (independent hits) within a significant eQTL. Consequently, the eQTL analysis was repeated for each significant eQTL gene after additionally adjusting the linear regression model for the most significant variant identified for the eQTL gene. The procedure was reiterated until no additional significant variants were identified. In this analysis, a variant was regarded a significant independent eQTL for a given gene if the P-value associated with the regression slope was lower than 1 × 10 −6 . With this approach, we detected an additional 101 independent eQTL variants in 93 out of 1,959 liver eQTL genes ( Fig. 1, Supplementary Tables S2 and S3). Of note, our analysis could not replicate the AMD associated eQTL rs797037040 previously shown to influence the expression of TNFRSF10A in blood 14 . This is owed to the fact that neither this variant nor any variant in linkage disequilibrium (R > 0.4) to rs797037040 could be reliably imputed into the dataset.

Characterization of eQTL-variants.
We further localized all independent eQTL hits with regard to the transcription start site (TSS) of the affected gene (Fig. 2). We observed that the most significant eQTL variants were close to a respective TSS. Overall, 1,599 out of 2,060 (1,959 + 101) independent eQTL variants were within 100,000 base pairs of a nearest TSS, well in agreement with other studies 16,[25][26][27] .
We then evaluated the RegulomeDB 28 scores of eQTL variants ( Fig. 3A and Supplementary Table S4). As expected, eQTL variants (N = 183,872) were enriched in RegulomeDB classes one to four (P-values < 6.82 × 10 −09 ), which represent variants with likely regulatory properties while categories 5 and higher show minimal to no functional relevance. In addition, eQTL variants with the smallest P-values and additional secondary signals (independent hits, N = 2,040) revealed an even stronger enrichment in classes one to four compared to controls and compared to all eQTL variants (P-values from 1.72 × 10 −04 to 8.27 × 10 −11 ).
To further characterise each eQTL signal for its most severe functional consequence relative to a known gene structure, we applied Ensembl VEP 29,30 (Fig. 3B, Supplementary Table S5). Control variants were predominantly located upstream (49.22%) and downstream (49.09%) of known gene structures. Another 1.63% of the control variants were found in introns of genes. Less than 0.1% of the control variants were assigned to functional Study Schadt et al. 18 Schroeder et al. 19 Innocenti et al. 17 GTEx Start/Mid a Meta-analysis Mega-analysis  Table 1. Study and sample summary QC = quality control; a Omni 2.5 M for the first data release (GTEx start) and Omni 5 M for the mid-point release (GTEx mid). b After quality control the genotype files of the four studies were merged into a single file and variants, which did not overlap in-between datasets, were assigned missing. We only kept variants which were genotyped in at least 100 samples.  Table S5). Additionally, other predicted categories like UTR or coding region variants occurred more often (P-values < 1.72 × 10 −07 ). Taken together, our findings indicate that significant liver eQTL variants are more often localized within known gene structures and are likely regulatory variants as they are found within regions of transcription factor binding and open chromatin. In addition, the most significant variants are also the most likely functional variant in each eQTL. This is supported by findings that the most significant eQTL variants (i) show an increased level of enrichment in all relevant RegulomeDB score categories compared to all eQTL variants and (ii) are enriched within known gene structures such as introns or coding exons.
Liver eQTL in AMD. Finally, we investigated whether any of the 52 independent AMD associated variants reported by Fritsche et al. 3 coincides with the established liver eQTL. Out of 52 independent tag variants, only 31 variants had an allele frequency >5% and could be reliably imputed into our dataset. Remarkably, 8 of these 31 variants significantly affect 15 unique eQTL-genes (Q-Value < 0.05, Table 2).
Within the complement factor H (CFH) locus, several AMD associated variants appear to influence expression of CFH and CFH related genes (CFHR). Particularly, the independent hit variant rs10922109 (independent hit 1-1 in 3 ) tags a common deletion of CFHR1/CFHR3. Since the deletion of both genes is protective against AMD, the risk increasing allele results in elevated expression of the two genes (Table 2).
Notably, two genes involved in HDL metabolism, Cholesteryl ester transfer protein (CETP) and hepatic lipase (LIPC), were both significantly regulated by AMD associated variants (Table 2). Specifically, rs17231506 is highly correlated to rs3764261 (R² > 0.99), a variant that results in markedly increased HDL levels in blood 31 . According to our eQTL data, rs17231506 reduces the expression of CETP, in line with the observation that CETP deficiency or pharmacological inhibition leads to elevated serum HDL. Further, our eQTL data showed that rs2070895 (−250 G > A) increases the expression of LIPC and would be expected to be associated with decreased HDL blood 32 .
Finally, we identified additional AMD associated variants that potentially act as eQTL in liver. The AMD risk increasing allele of rs7803454 increases the expression of the paired immunoglobin like type 2 receptor alpha (PILRA) and beta (PILRB) genes. The resulting proteins are known to function as antagonists within the Tyrosine-protein phosphatase non-receptor type 6 (PTPN6) pathway 33 and have been implicated in both, AMD and Alzheimer's disease risk 34 . Interestingly, we did not detect any eQTL within the strongest AMD associated locus located on chromosome 10q26 (ARMS2/HTRA1).

Discussion
In this study, we have combined the genotypes and expression data of four previously published independent studies to further our understanding of the regulatory networks in liver tissue. Each individual study intended to identify new liver specific eQTL in order to elucidate the contribution of regulatory mechanisms on different diseases or traits. For example, Schadt et al. 18 were the first to explore eQTL in liver tissue and correlated their results to genome-wide association studies of seven different diseases. AMD was not among them. Innocenti et al. 17 and Schroeder et al. 19 followed a similar approach but concentrated on the reproducibility of eQTL, while the latter group additionally focused on genes involved in drug response pathways. GTEx analysed eQTL in 44 human tissues and aimed to explore the interplay of gene regulation across tissues. By merging these resources this is to our knowledge the largest study on liver eQTL to date and promises to provide novel insight into the role of genetic variation on gene expression in liver tissue. Combining several studies while jointly analysing the data has drastically increased the power to detect novel eQTL across the genome. The replication rates of eQTL detected in individual studies can be as low as 38.5% (Table 1), even with a stringent FDR threshold of 0.1%. An approach known as mega-analysis has further improved the power of our study to detect novel eQTL. This also revealed a higher replication rate of eQTL identified by individual studies. Although the gain in power attributable to a mega-analysis can depend on the type of study 23 , the mega-analysis approach allowed us to identify additional, independent signals in 5% of the significant eQTL.
Mapping identified eQTL-variants against known gene structures such as introns, coding or non-coding exons revealed that a large proportion of the identified eQTL variants is highly enriched in intronic and coding regions of genes, in line with previous results 13,16 , although such an enrichment may be specific for certain tissues 35 . Similarly, we have observed a strong enrichment of eQTL variants in RegulomeDB classes one to four representing known eQTL and expected regulatory variants. Since many eQTL are shared between tissues 20 , an enrichment in RegulomeDB class 1 (representing known eQTL) is not surprising. Nevertheless, we also observe a strong enrichment of eQTL variants in RegulomeDB classes two to four, representing variants in experimentally determined regulatory epigenetic elements. Importantly, hypothetic regulatory variants in RegulomeDB class 5 (characterized by either transcription factor binding or a peak of DNase hypersensitivity) are not enriched in the identified liver eQTL variants, greatly increasing confidence in the robustness of our results. Alternatively, variants in RegulomeDB class 5 could be variants with weaker regulatory effects and thus, our study might be underpowered to identify significant eQTL variants that are characterized by mapping to a weak epigenetic mark.
Strikingly, the observed enrichment in gene structures were more pronounced in the independent hits which represented the most significantly associated variants and, in addition, the most significantly associated secondary signals. This strengthen the notion that the variant showing the smallest P-value of association or correlation in a locus is a priori the most likely one to be the true causative mutation 36 . Alternatively, it is also possible that the functional allele of the variant with the smallest P-value is rather tagging several haplotypes that affect gene expression in the same orientation 37 . Therefore, in case a defined eQTL is of major interest, such a locus has to be dissected further by statistical means to identify all independent haplotypes carrying functional alleles 10 .
While the central nervous system and the retina are expressing complement genes, the liver is nevertheless the primary site of synthesis for circulatory complement proteins 38 . In addition, the liver plays a key role in lipid metabolism 39 , besides the complement cascade another pathway implicated in AMD pathology by epidemiological and genetic studies. We therefore investigated whether any of the top hits of a recent GWAS for AMD 3 are regulatory variants influencing gene expression in liver.
One of the most significant association signals for AMD resides within the CFH locus on chromosome 1 and represents a compound signal of two protective haplotypes tagged by the protective allele of the top variant 37 . One protective haplotype harbors a common deletion of the CFH-related genes 1 and 3 (CFHR1/3) 40 . The heterozygous deletion of both genes results in reduced levels of CFHR1/3 proteins in serum, while a homozygous deletion results in a complete absence of CFHR1/3 41,42 . In line with this, we found that the AMD risk increasing allele of rs6677604 is correlated to increased expression of both genes while the protective allele of rs6677604 (in strong linkage disequilibrium with the CFHR1/3 deletion) is correlated with decreased expression. In addition, the protective allele reduces the expression of other CFHR genes as well as the expression of the CFH gene. Since CFH and CFH-related genes share high sequence identity with each other, the expression values of the individual gene may not be distinguishable from the related gene by currently used high-throughput methods [43][44][45] . Indeed, we found that the gene expression values of CFH and CFH-related genes (CFHR1-5) are correlated in liver samples (R² between 0.1 and 0.5).
One important result of our study reveals that two AMD-associated signals near LIPC and CETP are significant eQTL, strongly implicating HDL metabolism and serum lipid levels in AMD pathogenesis. We observed that the AMD risk increasing allele of rs17231506 reduces CETP expression, likely resulting in elevated HDL levels in serum 46 . This is in line with the observation that HDL levels are elevated in AMD patients compared to controls [7][8][9] . Further, the risk increasing allele of rs2070895 near LIPC results in increased expression of LIPC, which is generally associated with reduced serum HDL levels 47 . A study by Burgess and Smith 48 also observed an AMD associated variant next to LIPC (rs261342) to be associated with decreased HDL serum levels 48 . This variant is in high linkage disequilibrium with rs2070895 (R 2 = 0.84) which was shown in our study to cause elevated LIPC expression in liver. Burgess   key regulators of HDL remodelling which might be essential for efficient delivery of lipids (e.g. fatty acids, carotenoids) into the retina and efflux of excess lipids out of the retina. Importantly, CETP and LIPC variants have been shown to have additive effects on cardiovascular risk with low CETP activity variants combined with low LIPC activity variants increased the risk 49 . Cardiovascular risk could therefore add additional pressure to select specific variant gene combinations in the aged AMD population that were protected from cardiovascular death. A similar line of thought emerged from another recent study, which found that a genetic score based on genome-wide significant variants for elevated HDL serum levels was higher in AMD patients, strongly suggesting that AMD patients have more alleles that increase HDL than controls 50 , in line with other studies 51,52 . Other confounding variables such as exercise, drugs or alcohol consumption or the occurrence of AMD in study participants are potentially influencing our eQTL analysis. However, the individuals in the study were largely below 60 years of age (404 out of 588) and thus AMD associated impairment such as an overly sedentary life style should play a minor role in confounding our analysis. Furthermore, this study included a diverse and large set of individuals across multiple studies, which should reduce the effect of confounding environmental factors, especially since AMD associated factors are not likely to significantly influence confounders such as alcohol consumption 53,54 or treatment with different, liver-metabolized drugs.

Conclusions
We present the currently most comprehensive eQTL analysis for liver tissue and report that 1,959 out of 24,123 investigated genes have at least one significant eQTL in liver. Significant eQTL variants are more frequently found within gene boundaries and are more enriched in RegulomeDB classes representing likely regulatory variants. Several of these liver eQTL overlap with genetic variants strongly associated with AMD at genome-wide significance. These findings underscore the validity of the eQTL approach to identify disease-associated functional variants and provide further confirmation that HDL metabolism is strongly involved in AMD aetiology. Nevertheless, it should be emphasized that further replication of our results in disease relevant tissues such as retina or RPE or other functional validation studies are warranted. Specifically, this could further validate our notion that HDL metabolism is, in addition to the complement cascade, a major pathway in AMD disease development.

Methods
Genotype data. The genotypes of the four studies were retrieved from the respective databases (Table 1).
Genotype quality control was performed for each study separately and, in addition, jointly after imputation. Since some studies reported only the zygosity of their samples at each variant (e.g. homozygosity: AA or BB; heterozygosity: AB), we first matched the reported alleles of each variant to the respective allele in the 1000 Genomes reference dataset to the Biomart 30,55 online database (http://grch37.ensembl.org/biomart/). Multi-allelic variants were excluded to avoid potential ambiguity. Next, for each study we extracted the genotypes of all samples at 30,000 randomly chosen variants from all autosomes. We also included the genotypes of all samples from the 1000 Genomes Project Phase 3 (release 20130502) 56 at the same variants and performed a PCA with the snpgd-sPCA function of the SNPRelate 57 package in R 58 . Since the haplotype structure can greatly vary between populations, we only included individuals clustering next to the European (EUR) reference individuals in the eQTL analyses ( Supplementary Fig. S1). We then compared the reference allele in the datasets to the reference allele in the European 1000 Genomes samples. Alleles were flipped when given on the opposite strand. We excluded variants whose reference allele frequency differed by more than 10% from the reference allele frequency of the 1000 Genomes European samples. Furthermore, we excluded variants that were (1) not on autosomes, (2) had a minor allele frequency of MAF < 0.05 or deviated significantly from Hardy-Weinberg equilibrium 59 (HWE, P < 1 × 10 −6 ) after applying the respective function in the VCFtools 60 . The individual genotype data sets were merged into a single VCF file. Variants which were not present in an individual study or were not genotyped in at least 100 samples were assigned missing in the respective individuals. Phasing and imputation was performed on the merged data, as accuracy of both algorithms increases with increasing sample sizes 61 . Phasing was performed with SHAPEIT2 and standard settings by supplying the imputed genotypes from the 1000 Genomes Phase 3 reference panel 62 . The same reference panel was used to conduct a whole genome imputation with IMPUTE2 63 at standard settings. Next, VCFtools was used to remove variants with a minor allele frequency < 5% and variants which showed evidence for a significant deviation from Hardy-Weinberg equilibrium (P < 1 × 10 −6 ). In addition, variants with an IMPUTE2 info score smaller than 0.4 considered to be of low quality 64 , were removed. Finally, the reference allele frequency of each study was compared against the reference allele frequency of all other studies ( Supplementary Fig. S2). Variants whose reference allele frequency differed by more than 15% between studies were excluded.
Specifics for each data set were as follows: The GTEx data were retrieved through dbGAP 65 (https://www.ncbi.nlm.nih.gov/gap, accession: phs000424. v6.p1). The positions of the variants were already reported based on the final hg19 build and thus, no additional lift-over was required.
Innocenti et al. 17 genotype information was retrieved from the GEO database 66 (accession code: GSE26105). The genotyping had been performed by the authors on an Illumina 610 Quad chip and the genotypes were encoded by each individual's zygosity status (homozygosity: AA, BB; or heterozygosity: AB). The hg19 coordinates as well as the respective alleles of the variants were retrieved from Ensemble by querying the Biomart online database with the respective dbSNP identifier.
The genotype information from Schroeder et al. 19 was retrieved from the GEO database (accession: GSE39036). The samples had been genotyped by the authors on an Illumina HumanHap300 chip and the genotypes were also encoded according to the individual zygosity status. The hg19 coordinates and alleles were retrieved from Ensemble as specified above. The genotypes from the Schadt et al. 18 study were retrieved from the Synapse database (accession: syn89614). The samples had been genotyped on either the Affymetrix 500k or the Illumina 650 Y genotyping chip. The genotype file included hg17 positions of each variant, a unique dbSNP identifier and both alleles of each individual. We initially removed variants without dbSNP identifiers and then used the program liftover 67 from the UCSC Genome Browser (https://genome.ucsc.edu/util.html) to retrieve the hg19 coordinates of each variant.
Gene expression data. The present study included the gene expression data from four independent studies. Three studies profiled gene expression by employing microarray platforms (Table 1) while one study used high-throughput transcriptome sequencing (RNA-Seq) for data generation. First, we remapped array probes to an in silico mRNA reference database based on Ensemble gene annotation 30 with the help of the ReAnnotator pipeline 68 . Only exome-matching probes showing less than five mismatches were retained in the data set. Probes mapping to multiple genes or overlapping with common variants (according to dbSNP release 142) were removed from the analysis 69 . Probes which measured the gene expression of the same gene, were merged by calculating the mean of all probes within a gene, weighted by the variance of the respective probe over all samples. Hence, probes with a higher variance contributed more to the overall transcript levels than probes with little variation across samples.
For each data set, we performed basic expression normalization and quality control. Briefly, the available expression values were log2-transformed and a PCA was performed with the prcomp function in R to detect potential outlier samples within the dataset. We merged replicate samples by taking the mean of all replicate values.
The expression data of the four studies were merged and missing expression values were imputed using the K-Nearest-Neighbour 70 method provided by the impute.knn function of the impute Bioconductor package 71 in R. Genes that were included in one study but could not be imputed into the other studies were removed. Differences between all individuals were evaluated by conducting a PCA on the gene expression data ( Supplementary  Fig. S3A-C). In addition, the expression values for each individual were plotted as a boxplot ( Supplementary  Fig. S3D-F). Due to substantial differences between datasets, we applied further normalisation steps. Initially, we performed a quantile normalisation with the normalize.quantiles function of the R package preprocessCore 72,73 . Since quantile normalization alone was not sufficient to normalize all studies, we adopted an empirical batch correction method called ComBat with the combat function from the sva package in R 74 . By supplying known batch effects to the function (i.e. the study labels), ComBat standardises the data gene-wise and then applies an empirical batch effect correction ( Supplementary Fig. S3C and F). The batch corrected expression values were used for the eQTL analyses, as no obvious bias of the single studies was noticeable.
Methods specific to the individual studies were as follows: Firstly, for the GTEx data expression values (release GTex-V6p) were downloaded from the GTEx Portal (http://www.gtexportal.org/home/). The levels of transcript expression were encoded as "reads per kilobase of transcript per million mapped reads" (RPKM). We added 0.001 to all RPKM values to perform a log 2 transformation of the data.
Secondly, the expression data from Innocenti et al. 17 were retrieved from the gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/, accession: GSE25935). The expression values were already background subtracted and transformed to the log 2 scale.
Thirdly, Schadt et al. (2008) provided a curated version of their data in the Synapse database (https://www. synapse.org/, accession: syn89614). As this study used an Agilent Custom 44k array, probe sequences were not openly available. In addition, not all samples had values for both genotype and gene expression data. The authors supplied an annotation file which links probe IDs to Ensemble and RefSeq 75 identifiers. Expression values of probes were only used, if they were unanimously linked to a single Ensemble or RefSeq identifier. Furthermore, RefSeq identifiers were converted to Ensemble gene identifiers with the help of the Ensembl biomart tool 55 . A Shapiro-Wilk test 76 revealed that raw values larger than 2 or smaller than −2 values are likely outliers. Thus, all of these were set to missing.
Finally, expression values from Schroeder et al. 19 were retrieved from the GEO database (accession: GSE32504) as quantile normalized data. To retrieve probe sequences of the Illumina Human WG-6v2.0 chip for probe remapping, the illuminaHumanv2.db R package 77 was used. eQTL analysis. Linear regression analysis between gene expression values and imputed allele dosages was performed with Matrix eQTL 78 . Age, gender and the first five principal components of the genotype PCA were included in the models as covariates. We exclusively calculated local eQTL (variant-gene distance less than one million base pairs) due to limited power to perform distant eQTL analyses 15 .
Two approaches were adopted to jointly analyse eQTL. First, a classic meta-analysis was applied to the individual study results. The effect size (slope) and standard error of the effect size were estimated with Matrix eQTL for each study separately. Further, a random effects model implemented in the function MiMa 79 was applied to estimate the joint effect sizes and standard errors as well as the joint P-Values. The latter approach (mega-analysis) estimated local eQTL from the merged genotype and expression data directly. This approach also allowed us to search for novel independent eQTL for a gene by adjusting the linear regression model for the most significant eQTL variant for this gene. To account for multiple testing, the false discovery rate (FDR) was controlled to be smaller than 0.001. Thus, joint Q-Values were considered to be smaller than 0.001 for statistical significance. The score denotes the confidence that a certain variant is important for transcription factor binding or chromatin accessibility and thus gene regulation. Variants in classes one to four are deemed very likely regulatory variants, while variants in classes five to seven are less likely to influence gene expression. In addition, the Ensembl Variant Effect Predictor (VEP, www.ensembl.org/vep) was used to assign each eQTL variant to a functional consequence relative to known gene structures. The program predicted the most severe consequence per gene within a range of 1 Mbp up and downstream of each variant. For eQTL variants, only predicted consequences affecting the associated eQTL gene were evaluated. For the control variants, a single random consequence for a nearby gene was chosen.
Ethics approval and consent to participate. This study used data of four public datasets. For further specifics on the respective ethics approvals, we refer to the single study publications.
Data availability statement. All data are available in public databases as detailed in the methods section.