Expression Quantitative Trait Loci Acting Across Multiple Tissues Are Enriched in Inherited Risk for Coronary Artery Disease

In the STAGE study, 7 vascular and metabolic tissues of well-characterized CAD patients were sampled during coronary artery bypass grafting as described.¹⁰ The Ethical committee of the Karolinska Hospital approved the study, and patients gave written consent (Dnr 004-02). Patients were included if they were eligible for coronary artery bypass grafting and had no other severe systemic diseases (eg, widespread cancer and active systemic inflammatory disease).

STAGE tissue samples were obtained from atherosclerotic arterial wall (AAW), internal mammary artery (IMA), liver, skeletal muscle (SM), subcutaneous fat (SF), and visceral fat (VF)¹⁰ and stored at –80°C in RNAlater (Ambion). Preoperatively, patients were examined and fasting whole blood (WB) was obtained for DNA and RNA isolations and biochemistry. DNA and RNA were isolated with Qiagen Mini/Midi-Kits (Ambion). DNA from 109 patients (Table 1) with sufficient quantities and qualities (≥1 μg, 1.7>260/280>1.9 with Nanodrop, Agilent) was genotyped with GenomeWideSNP_6 arrays (Affymetrix). Allele frequencies for 909 622 SNPs were determined with the Birdseed algorithm in Affymetrix Power Tools (v 1.14.2); 530 222 autosomal SNPs with call rates of 100% and minor allele frequency >5% and in Hardy-Weinberg equilibrium (P>10^–6; quality controlled [QC] SNP set) were used for downstream analysis. Custom-made HuRSTA-2a520709 arrays (Affymetrix) were used for gene expression profiling of total RNA samples (≥5 μg, 1.95>260/280>2.05 with Nanodrop, Agilent) from the 109 genotyped patients (AAW, n=68; IMA, n=79; liver, n=77; SM, n=78; SF, n=63; VF, n=88, and WB, n=102). Robust multiarray average was used for background subtraction, normalization, and summarizing of raw microarray data (Affymetrix Power Tools, v 1.14.2). A custom-made chip description file was used to match 381 707 probes on the array to 19 610 probe sets for unique genes (to avoid cross-hybridization between alternative transcripts) according to the hg19 human genome assembly.

Each expression trait (probe set signal) was tested for association with the allele frequency of each genotyped SNP using the Kruskal–Wallis test and false-discovery rates (FDRs) to control family-wise error rates.¹³ First, an empirical null distribution was estimated for each tissue by calculating association P values between SNPs and expression traits in data sets with randomized sample labels (patient identification numbers). Next, the FDR at level α was controlled by determining the P value threshold P_α such that the ratio of the average number of associations with P value <P_α under the null distribution to the number of associations with P value <P_α in the real data was equal to α. We found that the empirical null distribution was stable and 3 permutation runs per tissue were sufficient for convergence, indicating the distributions among the expression traits were largely exchangeable. At an FDR of 15%, the nominal genome-wide P-value thresholds for the Kruskal–Wallis test were 2.9×10⁻⁷ (AAW), 2.9×10⁻⁷ (IMA), 6.0×10⁻⁷ (liver), 3.5×10⁻⁷ (SM), 4.0×10⁻⁷ (SF), 2.3×10⁻⁷ (VF), and 6.0×10⁻⁷ (WB). Cis-acting expression single-nucleotide polymorphisms (eSNPs) were defined as SNPs within 1 Mb of the transcription start or end sites of the affected gene. All other eSNPs were defined as trans-acting. Because multiple testing corrections considering all SNPs (≈10¹⁰ pairwise tests) are too conservative for identifying cis-eSNPs (≈7×10⁶ pairwise tests), the above procedure was repeated, this time testing only SNPs within 1 Mb of the transcription start or end site for each gene. At the threshold of FDR=15%, nominal P-value thresholds for the Kruskal–Wallis test for cis-acting eSNPs were 2.5×10^–4 (AAW), 2.9×10^–4 (IMA), 5.5×10^–4 (liver), 3.0×10^–4 (SM), 3.0×10^–4 (SF), 2.4×10^–4 (VF), and 6.7×10^–4 (WB). This resulted in a total of 29 530 cis-eSNPs for 6450 genes, corresponding to 8156 eQTLs (the eSNP most strongly associated with each gene), and 1494 trans-acting eSNPs for 1222 genes (Table 1; Methods in the Data Supplement). Of 29 530 cis-eSNPs, 15 012 eSNPs were identified in WB, which, at least in part, likely is attributed to the relatively larger sample number of this tissue (n=102). However, on the whole, variation in STAGE sample sizes (n=63–88 besides WB) seemed to have subtle, if any, effect on the number of identified eSNPs. For regional plots, LocusZoom¹⁴ was used for selected CAD/MI risk loci, plotting for each SNP in the QC set its lowest association P value over all genes in all tissues and using the Europe population, hg19 genome assembly, and 1000 Genomes¹⁵ to calculate linkage disequilibrium (LD).

We analyzed 555 mRNA profiles (considered as 19 610-dimensional vectors) by multidimensional scaling using the Euclidean distance and Sammon’s nonlinear mapping criterion. Tissues were also clustered according to eSNP overlap by constructing a phylogenetic tree from 29 530-dimensional vectors recording the presence/absence of an eSNP in a tissue, using the neighbor-joining method with average linkage and the Jaccard distance measure.

Enrichments of eSNPs at genomic regions with certain features were compared to the bulk (n=963) of verified risk loci¹⁶ identified in GWA studies³ and present in our QC SNP set, using chromosome-bound circular permutations.¹⁷ The null hypothesis of the permutation method is that eSNPs and non-eSNPs are equally likely to be located (ie, found) within each genomic annotation. The real data were compared against 20 000 permuted genomes generated from the QC SNP set, which preserved the internal structure of the genome in terms of relative distance between SNPs and the observed clustering of genomic features while analyzing the LD structure (LD partners, r²>0.9) around all SNPs (Methods in the Data Supplement).¹⁷ The permuted genomes were used to calculate the odds ratios of risk enrichment/depletion of the real data in genic and regulatory regions, evolutionarily conserved regions, and regions with distinctly different histone modification patterns (Table I in the Data Supplement). Human umbilical vein endothelial cells were the most CAD-relevant cell type for which histone modification data were available.

To examine eQTLs and tissue-defined subsets thereof for relative enrichment in CAD/MI risk,^18,19 we used the MIGen¹¹ and WTCCC¹² GWA data sets. The queried eQTL set (real set >10 eQTLs; Table II in the Data Supplement) was matched against SNPs of the GWA data set and expanded with 1000 Genomes¹⁵ to include SNPs in that data set in strong LD (r²>0.8) with eQTLs. The expanded set was compared with 10 000 SNP sets of equal size randomly sampled from each GWA data set with matching chromosomal distribution. The fold enrichment was computed from the percentage of disease-associated SNPs (P<0.05) in the expanded set compared with the average percentage of disease-associated SNPs in the random sets. The enrichment P value was determined from the percentage of random sets with the same fold enrichment as the real set (or higher; Methods in the Data Supplement).

For highly risk-enriched sets of eQTLs, many associated gene sets was identified. Such gene set was identified by first seeking genes associated in cis or in trans with each of the eQTLs of interest across the 7 tissues at a nominal P value for the Kruskal–Wallis test (P<1×10^–3). Next, highly coexpressed genes active in the same tissue as eQTL-gene pair was identified forming the final gene set (Pearson correlation coefficient, >0.85). The threshold parameters were those yielding optimal gene set sizes from all combinations of association P-value thresholds of 5×10^–4, 1×10^–3, and 5×10^–3 and Pearson correlation threshold of 0.8, 0.85, and 0.9. eQTLs and associated gene sets in a tissue were assessed for molecular functions in Gene Ontology categories.²⁰ Functional-enrichment P values were determined with Fisher exact test and Benjamini and Hochberg FDR corrections for multiple testing.

As previously described,²¹ THP-1 monocytic cells were differentiated into macrophages using phorbol myristate acetate and thereafter transfected with small interfering RNA targeting each master regulator (one at a time) or mock-transfected (Table III in the Data Supplement). Next, the cells were incubated with acetylated low-density lipoproteins (50 μg/mL) for 48 hours allowing foam cells to form. The cholesterol-ester accumulation in the foam cells were determined by measuring and subtracting free cholesterol from total cholesterol adjusted for protein concentration (Methods in the Data Supplement).

The STAGE cohort consisted of typical CAD cases with well-managed hyperlipidemia and hypertension (>90% on β-blocker and plasma-lipid-lowering regimens); 20% had diabetes mellitus (Table IV in the Data Supplement). The expression patterns of the 19 610 genes across 7 tissues in the STAGE study segregated the tissues into well-defined groups, reflecting unique tissue functions (Figure 1A). Gene expression patterns in AAW and nonatherosclerotic IMA were closely related but, importantly, well segregated; gene expression patterns in SF and VF were more intermixed (Figure 1A; Figures I and II in the Data Supplement).

The allelic distributions of 29 530 SNPs within 1 Mb of the transcription start site (TSS) were associated with mRNA levels of 6450 unique genes (cis-eSNPs [eSNPs]; FDR 15%; Methods in the Data Supplement; Table 1; and Figure III in the Data Supplement). Of the identified 29 530 eSNPs, 22 101 were identified in only one tissue (tissue-specific eSNPs) and 7429 in ≥2 tissues (multitissue eSNPs). Similar proportions of tissue-specific versus shared eSNPs have been reported.^22–25 Also identified were 1494 trans-eSNPs (>1 Mb from the TSS), mainly in WB (38%). Within these sets, we identified 39 eSNPs matching SNPs previously associated with CAD/MI risk in GWA studies,⁴ whereof 21 exact matches and 18 in strong LD (r²>0.8; Table V and Figure IV in the Data Supplement). On average, the multitissue eSNPs were closer to the TSS than the tissue-specific eSNPs (Figure 1D and 1E), as reported.^22,26 Whether the multitissue eSNPs were present in 2 or more tissues did not affect their positions relative to the TSS (Figure 1D). The notion that eSNPs more distal to the TSS have tissue-specific effects on gene expression is consistent with the fact that tissue-specific enhancer regions are commonly more distal to the TSS. Similarly, the notion that eSNPs close to the TSS have tissue-independent effects on gene expression (across 2 or more tissues) is consistent with the idea that the TSS of a gene is similar across tissues.

eSNPs in arterial wall samples (AAW, IMA) and fat deposits (VF and SF) formed a common cluster, as did those in liver and WB; eSNPs in SM formed an isolated cluster (Figure 1B). In fact, ≤60% of the eSNPs in AAW, IMA, VF, and SF were shared (Figure 1C). Thus, in CAD, diseased and healthy arteries seem to share most regulatory SNPs with fat depots, fewer with liver and WB, and fewest with SM. Venn diagrams for eQTLs in all combinations of 4 tissues are shown in Figure V in the Data Supplement.

The strongest eSNP association of a gene is referred to as its eQTL. To assess the importance of different sets of tissue-specific and multitissue eQTLs for inherited risk of CAD/MI, we used the MIGen¹¹ and WTCCC¹² GWA data sets to analyze risk enrichment (Table II in the Data Supplement). According to both data sets, the entire set of cis-eQTLs was risk enriched by 1.15- to 1.2-fold (Table II in the Data Supplement). The set of trans-eQTLs was not risk enriched (not shown).

CAD is a complex disease that is particularly dependent on the molecular status of the metabolic tissues that drive atherosclerosis and thrombosis in the vascular bed. Thus, eQTLs present only in one tissue (tissue-specific) may have different roles from multitissue eQTLs, affecting gene expression across several tissues in mediating CAD risk. We examined this possibility by assessing risk enrichments of different sets of eQTLs present in single or multiple tissues (>2 tissues). According to MIGen GWA data set, only tissue-specific eQTLs in SM, SF, and WB carried CAD/MI risk (Figure 2A, red bars). However, for multitissue eQTLs, all tissues but SM and VF were risk enriched (Figure 2B, red bars). According to the WTCCC, the risk-enrichment levels for both tissue-specific and multitissue eQTLs were higher, but multitissue eQTLs were again found to be more risk enriched than tissue-specific eQTLs (Figure 2A and 2B, blue bars).

Because multitissue eQTLs (>2 tissues) were more risk enriched than tissue-specific eQTLs, we asked whether the CAD/MI risk of eQTLs increases with the number of tissues in which these eQTLs affect gene expression (from 2 to 7 tissues). Indeed, this proved to be the case: risk enrichment peaked in combinations of 4 to 6 tissues (1.2-fold, P<8.0×10⁻⁴) according to the MIGen data set and in combinations of 6 tissues (2.7-fold, P<2.0×10⁻⁸⁶) according to the WTCCC data set (Figure 2C). Thus, eQTLs that affect gene expression across 4 to 6 tissues seem to be the most relevant as carriers of CAD risk.

When we instead only considered the highest enrichment for a combination of specific tissues (from 2 to 7), the risk-enrichment scores according to both GWA data sets generally increased from tissue-specific eQTLs (Figure 2D, outer circle, ≤1.2-fold, P<4.0×10⁻⁹) to those active in combinations of 5 or more tissues (Figure 2D, inner circle, ≤8.5-fold, P<1.7×10⁻⁴³).

Most verified GWA hits are enriched at sites of regulatory elements,^17,27 suggesting that disease-risk loci affect disease by modifying regulatory hot spots rather than by modifying the expression of specific genes. Conceivably, eQTLs are enriched at regulatory sites. To investigate this possibility, we compared the regulatory element enrichment patterns of all the identified eSNPs (n=29 530) and the tissue-specific (n=22 101) and multitissue (n=7429) eSNPs with the regulatory element enrichment patterns of 963 established GWA loci¹⁶ (GWA-SNPs) in a range of regulatory element subcategories (Table I in the Data Supplement). Strikingly and as opposed to tissue-specific eSNPs (Figure VI in the Data Supplement), the regulatory element enrichment patterns for the multitissue eSNPs were significantly enriched, and more so than those for GWA-SNPs, regardless of the subcategory (genic and regulatory regions, evolutionary and conserved regions, or regions of histone modification in human umbilical vein endothelial cells; Figure 3A–3C). This finding strongly suggests that in contrast to tissue-specific eQTLs, and to a greater extent than GWA-SNPs, multitissue eSNPs/eQTLs affect gene expression across tissues mainly by modifying DNA regulatory elements at well-conserved genomic sites. This finding is consistent with their tissue-independent effects on gene expression (ie, being multitissue eQTLs) and possibly with their greater importance as carriers of CAD risk.

The finding that eQTLs shared across combinations of tissues (>4) were more enriched with inherited CAD/MI risk than tissue-specific eQTLs and those shared across fewer tissues (<4) suggests that there might be specific subsets of highly shared eQTLs that are particular important carriers of CAD risk. We therefore examined eQTLs active across specific combinations of the 7 tissues. Twelve eQTLs in one combination of 5 tissues (AAW, liver SF, VF, and WB) and 30 eQTLs in a combination of 6 tissues (IMA, liver, SM, SF, VF, and WB) were found highly enriched in CAD/MI risk according to the MIGen (8.5-fold, P=1.6×10⁻⁴³) and WTCCC (3.6-fold, P=2.2×10⁻³⁰) GWA data sets (Tables VI and VII in the Data Supplement). Combined, these 42 eQTLs were found to be risk enriched by ≤7.3-fold (7.3-fold, P=7.7×10⁻²⁰ in WTCCC and 2.3-fold, P=3.1×10⁻⁶ in MIGen). To validate the 42 eQTLs in a broader set of GWA studies of CAD, we interrogated the GWA meta-analysis of CAD undertaken by the Coronary Artery Disease Genome Wide Replication And Meta-analysis (CARDIoGRAM) consortium.⁴ The CARDIoGRAM data set consists of ≈22 000 CAD cases and ≈60 000 controls and includes the WTCCC and MIGen data sets, which contribute 14% of the cases and 7% of the controls.⁴ Using this GWA data set, we found that the 42 eQTLs were risk enriched by 4.2-fold (P<1.1×10⁻⁵⁴). The CAD associations for the individual SNPs in this set are shown in Table VIII in the Data Supplement.

To assess downstream effects on gene expression governed by the multitissue eQTLs, we first identified genes associated in cis or in trans with each of the 42 eQTLs (Figure VII in the Data Supplement; P<10^–3 by Kruskal–Wallis test) and then identified highly coexpressed genes downstream of the master regulators (absolute Pearson correlation coefficient, >0.85). In this fashion, 29 gene sets (with ≥30 genes in gene set in different tissues) for 16 of the 42 eQTLs, termed as master regulator eQTLs, involving 19 cis-regulated genes, termed as master regulator genes, were identified. According to the CARDIoGRAM GWA data set, these 16 eQTLs were further enriched with CAD risk (8.4-fold, P=5.5×10^–51). To put the fold enrichment of risk in perspective, we tested the risk enrichment for the 153 CAD risk SNPs established by GWAS² and found it to be 10-fold (P<10^–99) according to the CARDIoGRAM GWA data set.

To further evaluate the role of these 16 master regulatory eQTLs, we first used Gene Ontology²⁰ to assess the functional characteristics of the 29 gene sets. Nineteen gene sets returned significant functional enrichments that without exception revealed biological processes relevant to CAD: T-cell differentiation, regulation of cholesterol efflux, oxidation–reduction process, steroid metabolic process, blood vessel development, respiratory electron transport chain, and cell cycle (Table 2; Tables VI, VII, and IX in the Data Supplement).

We also examined whether clustering patients according to expression levels within each of the 29 gene sets would result in patient subgroups associated with CAD traits or clinical phenotypes. Remarkably, as many as 27 of the 29 gene sets clustered the STAGE patients according to ≥1 CAD trait or clinical phenotype (Table 2; Table IX in the Data Supplement). Heat maps of patient clusters within each of the 29 gene sets are shown in Figure VII in the Data Supplement.

In view of the master regulatory role of the 16 eQTLs in CAD, it is interesting to know more about their cis-regulated genes (master regulators). Three of the 4 master regulators with gene sets found to be associated with the primary CAD trait—the extent of coronary atherosclerosis assessed from preoperative angiograms in the STAGE patients (CAD score in Table 2)—were identified in abdominal fat (VF): G3BP1, PSORSIC3, and FLYWCH1. The gene sets of G3BP1 (a DNA-unwinding enzyme) and PSORSIC3 (a psoriasis-susceptibility gene) were highly enriched in cell-cycle genes (P<1×10⁻²¹ and P<1×10⁻⁴⁸). The gene set of FLYWCH1 (a zinc finger–containing protein) was instead enriched in oxidation–reduction processes (P<1×10⁻¹¹).

A blood (WB) gene set (n=53) with the master regulator SNAPIN (SNAP-associated protein) was also found to be associated with CAD score and enriched in genes active in T-cell differentiation (P<1×10⁻⁴; Table 2). This SNAPIN gene set was also associated with inherited diabetes mellitus, years of smoking, and physical activity in STAGE patients (Table 2; Table IX in the Data Supplement). Another SNAPIN gene set, identified in abdominal fat (VF), was found associated with low-density lipoprotein cholesterol levels (Table 2). For yet 2 other SNAPIN-associated gene sets (both identified in the arterial wall [IMA]) that were highly enriched in steroid metabolic processes (P<1×10⁻²⁰), one set was found associated with C-reactive protein levels of the STAGE patients and the other with the diagnosis of allergy (Table 2). Master regulator eQTL genes related to other phenotypes of the STAGE patients, including plasma levels of lipids, glucose, and C-reactive protein, are presented in Table IX in the Data Supplement.

The master regulatory genes G3BP1, PSORSIC3, FLYWCH1 (in VF), and SNAPIN (in WB) were all associated with CAD score in STAGE patients, suggesting these genes are important regulators of the formation of atherosclerosis. A key event of atherosclerosis is the formation of foam cells. To validate these genes, we examined the effect of small interfering RNA targeting of G3BP1, PSORSIC3, FLYWCH1, and SNAPIN on cholesterol-ester accumulation in THP1-foam cells. Besides SNAPIN, small interfering RNA targeting of G3BP1, PSORSIC3, and FLYWCH1 had pronounced effects on cholesterol-ester accumulation in THP-foam cells (Table 3; Table III in the Data Supplement).