Cardiovascular disease (CVD), which ultimately damages heart muscle, is a leading cause of death worldwide (WHO, 2018). CVD encompasses a range of pathologies including myocardial infarction (MI), where ischemia or a lack of oxygen delivery to energy-demanding cardiomyocytes results in cellular stress, irreparable damage, and cell death. Genome-wide association studies (GWAS) have identified hundreds of loci associated with coronary artery disease (Nikpay et al., 2015), MI, and heart failure (Shah et al., 2020), indicating the potential contribution of specific genetic variants to disease risk. Most disease-associated loci do not localize within coding regions of the genome, often making inference about the molecular mechanisms of disease challenging. That said, because most GWAS loci fall within non-coding regions, these variants are thought to have a role in regulating gene expression. One of the main goals of the Genotype-Tissue Expression (GTEx) project has been to bridge the gap between genotype and organismal level phenotypes by identifying associations between genetic variants and intermediate molecular level phenotypes such as gene expression levels (GTEx Consortium et al., 2017). The GTEx project has identified tens of thousands of expression quantitative trait loci (eQTLs); namely, variants that are associated with changes in gene expression levels, across dozens of tissues including ventricular and atrial samples from the heart. However, the eQTLs reported by GTEx explain a modest proportion of GWAS loci, and while increasing the diversity of tissues and sample sizes will enable further insight, orthogonal approaches also need to be considered.

It is becoming increasingly evident that many genetic variants that are not associated with gene expression levels at steady state, may be found to impact dynamic programs of gene expression in specific contexts. This includes specific developmental stages (Cuomo et al., 2020; Strober et al., 2019), or specific exposure to an environmental stimulus such as endoplasmic reticulum stress (Dombroski et al., 2010), hormone treatment (Maranville et al., 2011), radiation-induced cell death (Smirnov et al., 2012), vitamin D exposure (Kariuki et al., 2016), drug-induced cardiotoxicity (Knowles et al., 2018), and response to infection (Alasoo et al., 2018; Barreiro et al., 2012; Çalışkan et al., 2015; Kim-Hellmuth et al., 2017; Manry et al., 2017; Nédélec et al., 2016). The studies of context-specific dynamic eQTLs highlight the need to determine the effects of genetic variants in the relevant environment. Therefore, if we are to fully understand the effects of genetic variation on disease, we must assay disease-relevant cell types and disease-relevant perturbations. Most of the aforementioned studies were performed in whole blood or immune cells, which means that there are many cell types and disease-relevant states that have yet to be explored.

With advances in pluripotent stem cell technology, we can now generate otherwise largely inaccessible human cell types through directed differentiation of induced pluripotent stem cells (iPSCs) reprogrammed from easily accessible tissues such as fibroblasts or B-cells. One of the advantages of iPSC-derived cell types as a model system is that the environment can be controlled, and thus we can specifically test for genetic effects on molecular phenotypes in response to controlled perturbation. This is particularly useful for studies of complex diseases such as CVD, which result from a combination of both genetic and environmental factors.

The heart is a complex tissue consisting of multiple cell types, yet the bulk of the volume of the heart is comprised of cardiomyocytes (Donovan et al., 2019; Pinto et al., 2016), which are particularly susceptible to oxygen deprivation given their high metabolic activity. iPSC-derived cardiomyocytes (iPSC-CMs) have been shown to be a useful model for studying genetic effects on various cardiovascular traits and diseases, as well for studying gene regulation (Banovich et al., 2018; Benaglio et al., 2019; Brodehl et al., 2019; Burridge et al., 2016; de la Roche et al., 2019; Ma et al., 2018; McDermott-Roe et al., 2019; Panopoulos et al., 2017; Pavlovic et al., 2018; Ward and Gilad, 2019).

In humans, coronary artery disease can lead to MI (Dzau et al., 2006) which results in ischemia and a lack of oxygen delivery to energy-demanding cardiomyocytes. Given the inability of cardiomyocytes to regenerate, this cellular stress ultimately leads to tissue damage. Advances in treatment for MI, such as surgery to restore blood flow and oxygen to occluded arteries, have improved clinical outcomes. However, a rapid increase in oxygen levels post-MI can generate reactive oxygen species leading to ischemia-reperfusion (I/R) injury (Giordano, 2005). Both MI and I/R injury can thus ultimately influence the amount of damage in the heart. iPSC-CMs allow us to mimic the I/R injury process in vitro by manipulating the oxygen levels that cardiomyocytes are exposed to in vivo.

We thus designed a study aimed at developing an understanding of the genetic determinants of the response to a universal cellular stress, oxygen deprivation, in a disease-relevant cell type, mimicking a disease-relevant process. To do so, we established an in vitro model of oxygen deprivation (hypoxia) and re-oxygenation in a panel of iPSC-CMs from 15 genotyped individuals (Banovich et al., 2018). We collected data for three molecular level phenotypes: gene expression, chromatin accessibility, and DNA methylation to understand both the genetic and regulatory responses to this cellular stress. This framework allowed us to identify eQTLs that are not evident at steady state, and assess their association with complex traits and disease.

We differentiated iPSC-CMs from iPSCs of 15 Yoruba individuals that were part of the HapMap project (Banovich et al., 2018). To obtain a measure of variance associated with the differentiation process, and to more effectively account for batch effects, we replicated the iPSC-CM differentiation from three individuals three times, yielding 21 differentiation experiments in total. The proportion of iPSC-CMs in each cell culture was enriched by metabolic purification (see Materials and methods). On Day 20, cardiomyocyte differentiation cultures from each individual were split into sub-cultures for each of the four subsequent oxygen conditions and for assessment of cardiomyocyte purity (Figure 1—figure supplement 1). iPSC-CMs were matured by electrical pulsing and maintenance in cell culture for 30 days. On Day 30 (+/- 1 day), the median cardiomyocyte purity from one representative sub-culture across individuals was 81% (40–97% range), determined by flow cytometry as the proportion of cells that were positive for the cardiac-specific marker, TNNT2 (Figure 1—figure supplement 2; Supplementary file 1; see Materials and methods).

We studied the response of the iPSC-CMs to hypoxia and re-oxygenation (Figure 1A). To do so, we first cultured the iPSC-CMs at oxygen levels that are close to physiological oxygen levels (10% oxygen - Condition A) for 7 days. We then subjected the iPSC-CMs to 6 hours of hypoxia (1% oxygen - Condition B), followed by re-oxygenation for 6 hours (10% oxygen - Condition C), or 24 hours (10% oxygen - Condition D) as previously described (Ward and Gilad, 2019). Oxygen levels were reproducibly controlled in cell culture (Figure 1B, Supplementary file 1). In order to determine whether the cardiomyocytes were affected by the changes in oxygen levels, we measured the enzymatic activity of released lactate dehydrogenase throughout the experiment, as a proxy for cytotoxicity. We also measured released BNP, a clinical marker of heart failure (Maeda et al., 1998). As expected, both cytotoxicity (p=0.01, Figure 1—figure supplement 3A–B) and BNP (p=5×10⁻⁶, Figure 1—figure supplement 3C) levels increased following hypoxia and long-term re-oxygenation.

Figure 1 with 3 supplements see all

Download asset Open asset

With this system established, we sought to understand the contribution of the global gene regulatory response to the molecular and cellular response to hypoxia and re-oxygenation. To do so, we collected global gene expression data (using RNA-seq; n = 15 individuals), chromatin accessibility data (using ATAC-seq; n = 14), and DNA methylation data (using the EPIC arrays; n = 13; Figure 1C) in each condition. With these data we studied both the gene regulatory response to oxygen perturbation, as well as the interaction of the response with the underlying genotype of the assayed individuals.

Gene expression changes in response to hypoxia and re-oxygenation

We first sought to identify those genes important for regulating the response by analyzing the gene expression (RNA-seq) data. We processed samples in batches as described in Supplementary file 1 and mapped and filtered sequencing reads to prevent allelic mapping biases (Figure 2—figure supplement 1; Supplementary file 1; see Materials and methods; van de Geijn et al., 2015). We observed that one sample (18852A) was a clear outlier when comparing read counts for 18,226 autosomal genes across all samples, and thus excluded it from further analysis (Figure 2—figure supplement 2). We filtered out genes with low expression levels (see Materials and methods) to yield a final set with data from 12,347 expressed genes (see Materials and methods). We performed a number of correlation-based analyses using the data from the technical replicates (Figure 2—figure supplement 3), and confirmed that the quality of the data is high and that, in line with our flow cytometry data, our iPSC-CMs express a range of cardiomyocyte marker genes across conditions including MYH7 and TNNT2 (Figure 2—figure supplement 4).

We took advantage of the fact that we have replicate experiments from three individuals to correct the data for unwanted variation (see Materials and methods; Risso et al., 2014). Following this procedure, our samples clustered both by oxygen level and individual (Figure 2—figure supplement 5). To identify genes that respond to hypoxia and re-oxygenation, we first tested for differential expression between pairs of conditions using a linear model with a fixed effect for ‘condition’, a random effect for ‘individual’, and four unwanted factors of variation, learned from the data, as covariates. At an FDR of 10%, we identified thousands of genes that are differentially expressed between conditions (A vs. B = 4,983, B vs. C = 6,311; B vs. D = 6,792; A vs. D = 2,835; Figure 2 and Figure 2—figure supplement 6A). In order to identify a single set of genes which respond to hypoxia across conditions, we used Cormotif (Wei et al., 2015) to jointly model pairs of tests. This approach led to the classification of 2,113 genes (17% of all expressed genes) as responding to hypoxia, which we term ‘response genes’ (Figure 2, Figure 2—figure supplement 6B–C), and 9,949 genes that do not change their expression level throughout the experiment which we term ‘non-response genes’. Response genes are enriched for genes previously identified to respond to hypoxia in a Caucasian population of individuals (Chi-squared test; p<2.2×10⁻¹⁶; Ward and Gilad, 2019), and are highly enriched for gene ontologies in transcription-related processes (see Materials and methods, modified Fisher’s exact test; p=1×10⁻¹⁹).

Figure 2 with 6 supplements see all

Download asset Open asset

Dynamic eQTLs are revealed following hypoxia

Having established that oxygen stress initiates a transcriptional response affecting thousands of genes, we sought to identify eQTLs, either before or after oxygen stress. Using the combined haplotype test (CHT), an approach that leverages allele-specific information in small sample sizes (see Materials and methods; van de Geijn et al., 2015), we identified 1,886 SNPs which associate with gene expression (eSNPs) in at least one condition resulting in 1,573 genes with eQTLs (eGenes) in at least one condition (q-value <0.1; A: 613; B: 564; C: 564 and D: 464; Figure 3A–B; Supplementary file 2). Given that cardiomyocytes were split from a single differentiation culture for each condition, we do not expect cell composition to bias the identification of eQTLs in each condition. Indeed, the distribution of eQTL effect sizes is similar across conditions (Figure 3—figure supplement 1A), and the eQTL effect size values are correlated across conditions (Figure 3—figure supplement 1B). In addition, we included the same number of individuals in each condition, RNA extraction batches included all conditions from a given individual, RIN scores are similar across conditions, and the number of sequencing reads are similar across conditions (Supplementary file 1). Together, these results suggest that our overall power to detect eQTLs is similar across conditions. We refer to the 613 eGenes identified in condition A as ‘baseline eGenes’.

Figure 3 with 3 supplements see all

Download asset Open asset

Our goal was to identify dynamic eQTLs, which are either revealed or suppressed as the cells transition between conditions. Due to the small sample size of our study, we have incomplete power to detect eQTLs in any condition; thus, a naive comparison of eQTLs classified as ‘significant’ across conditions will result in an over-estimation of the number of dynamic eQTLs. Indeed, the p-values of genes whose expression levels are not significantly associated with a SNP in any condition deviate from the expected uniform distribution (Kolmogorov–Smirnov test; p<2.2×10⁻¹⁶; Figure 3—figure supplement 2A). To address this challenge, we first considered eQTLs identified using a q-value <0.1 in at least one condition, and visualized the p-value distributions of the corresponding eQTL associations in all other conditions. These p-value distributions are expected to be uniform if we had complete power to detect eQTLs in any condition (because in that theoretical case, even a naive comparison of eQTLs classified as ‘significant’ across conditions will result in the identification of true condition-specific eQTLs). Due to incomplete power, this is clearly not the case (KS test; p<2.2×10⁻¹⁶; Figure 3—figure supplement 2B); however, this distribution allowed us to choose a lenient secondary p-value-based cutoff, where values deviate from the uniform distribution, to classify dynamic QTLs (p=0.15; Figure 3—figure supplement 2B).

We specifically focused on two dynamic scenarios. First, we defined suppressed eQTLs, as eQTLs that are identified in condition A at a q-value <0.1 but not in any of the other conditions, with a p-value greater than 0.15 (37 instances; Figure 3C). Second, we defined induced eQTLs as eQTLs identified in conditions B, C, or D at a q-value <0.1, but not in A, with a p-value greater than 0.15 (330 instances; Figure 3C). This set of 367 dynamic eQTLs corresponds to 328 unique dynamic eGenes (see Materials and methods) and includes ZNF845 and RFC2 (Figure 3E) as well as PPARGC1, C8orf82 and CELF1 (Figure 3—figure supplement 3). While our choice of the particular statistical cutoffs is somewhat arbitrary, we can evaluate the false discovery rate associated with our chosen cutoff. Based on the p-value distributions of the corresponding eQTL associations in all other conditions, we estimate that our approach to classify dynamic eQTLs is associated with a false discovery rate of 48%. The relatively high FDR associated with our choice of statistical cutoffs does not indicate that these loci are not eQTLs; rather it means that if we had a larger sample size, roughly half of our dynamic eQTLs should have been classified as eQTLs in more conditions, potentially in all of them. As expected, using a more stringent p-value cutoff of 0.9, we find fewer dynamic eQTLs (19) and a lower false discovery rate of 7%.

We next wanted to determine whether the dynamic eQTLs we identified in iPSC-CMs are also eQTLs in primary heart tissue. To do so, we compared our 367 dynamic eSNPs and eSNPs regulating the same gene in all four conditions (‘shared eQTLs’, n = 20) to eQTLs identified in left ventricle heart tissue and 49 tissues assayed from hundreds of individuals in the GTEx study (GTEx Consortium et al., 2017). Thirty-two of 326 dynamic eSNP-eGene pairs tested in GTEx are eQTLs in heart tissue (9.8%), and 140 are eQTLs in at least one other assayed tissue (42.9%) demonstrating that many of these genes have context-dependent inter-individual variation in expression. Six of 19 shared eSNP-eGene pairs tested in GTEx are eQTLs in heart tissue (31.6%), and six are an eQTL in at least one other tissue (31.6%). Shared eSNPs are therefore more likely to be eSNPs in heart tissue than dynamic eSNPs (Chi-squared test; p=0.01). To further investigate eQTL concordance, we compared the eQTL effect size of our dynamic and shared eQTLs to the eQTL effect size in heart tissue determined by the GTEx consortium. Shared eQTLs overlapping heart tissue eQTLs (n = 6) have a higher concordance of effect than dynamic eQTLs overlapping heart tissue eQTLs (n = 32; Spearman correlation = 0.73 vs. 0.12) suggesting that our perturbation study has revealed novel eQTL effects.

Chromatin accessibility changes following hypoxia and re-oxygenation

We next asked whether the hundreds of transcription factor expression changes following hypoxic stress are accompanied by global chromatin accessibility changes. To examine this, we performed ATAC-seq experiments in all four conditions to identify regions of open chromatin (we were only able to collect these data from 14 of the 15 individuals; Figure 1—figure supplement 1, see Materials and methods). We filtered the ATAC-seq reads to include only those reads that map unambiguously to the nuclear genome (see Materials and methods, Figure 4—figure supplement 1). We identified a set of open chromatin regions in each sample, and merged samples across individuals within each condition. Genomic regions identified as accessible in each condition were then merged to yield a set of 128,672 open chromatin regions across conditions (with a median length of 312 bp). Regions with low read counts were filtered out, resulting in a final set of 110,128 regions. Analysis of various metrics revealed the data to be of good quality (Figure 4—figure supplements 2–3).

We sought to identify chromatin regions that are differentially accessible across pairs of conditions. Using a sensitive adaptive shrinkage-based approach with a False Sign Rate of 10% (Stephens, 2017), we could not detect changes in accessibility between baseline and hypoxia across 110,128 regions; however, we identified 831 differentially accessible regions (DARs) between hypoxia and short-term re-oxygenation (BC-DARs; 429 regions with increased accessibility and 402 with decreased accessibility), and 71 DARs between hypoxia and long-term re-oxygenation (BD-DARs; Figure 4A). Despite the fact that we do not identify any DARs between normoxia and hypoxia, there is a strong anti-correlation in the effect size between the normoxia (A) and hypoxia (B) conditions, and the hypoxia (B) and re-oxygenation (C) conditions across regions (Spearman correlation = −0.62; sign test p=4.6×10⁻¹⁴; Figure 4B) suggesting minor changes in accessibility in response to hypoxia that do not meet our significance threshold. Conversely, there is a strong correlation in effect sizes between hypoxia and short-term re-oxygenation (BC-DARs), and hypoxia and long-term re-oxygenation (BD-DARs; Spearman correlation = 0.74; Figure 4C), and 59 of the 71 BD-DARs are amongst the 831 BC-DARs (83%), suggesting that most regions have returned to baseline levels of accessibility by the first re-oxygenation condition. This includes a region within the intron of the FOXO1 gene, a master regulator of the oxidative stress response (Figure 4D). We therefore considered the 831 BC-DARs, henceforth DARs, in further analysis. These DARs represent 0.8% of the total number of accessible regions.

Figure 4 with 4 supplements see all

Download asset Open asset

Genomic features associated with differentially accessible regions (DARs)

We next wanted to determine what distinguishes DARs from constitutively accessible regions (CARs). To do so, we investigated three classes of genomic features: (1) promoter- and enhancer-associated marks, (2) transcription factor binding locations, and (3) underlying DNA sequence features. We found that DARs are more likely to overlap TSS than constitutively accessible regions (43% overlap vs. 11% overlap; chi-squared test, p<2×10⁻¹⁶; Figure 5A) suggesting that DARs may be involved in the gene regulatory response. Indeed, DARs are more likely to coincide with active histone marks in left ventricle heart tissue than constitutively accessible regions (H3K4me3: 57% overlap DARs vs. 24% overlap constitutively accessible regions; chi-squared test; p<2.2×10⁻¹⁶; H3K4me1: 86% overlap DARs vs. 52% overlap constitutively accessible regions; p<2.2×10⁻¹⁶; Figure 5B).

Figure 5

Download asset Open asset

To determine whether sequence-specific hypoxia-responsive transcription factors associate with differentially accessible chromatin, we integrated DARs with published chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) data for the well-studied hypoxia-inducible factors HIF1α and HIF2α (Schödel et al., 2011). A total of 234 of the 356 HIF1α ChIP-seq binding sites (66%), and 150 of the 301 HIF2α binding sites (50%) overlap with our set of 110,128 accessible chromatin regions. We found that these HIF1α and HIF2α binding sites are more likely to overlap the 831 differentially accessible regions than the 109,275 constitutively accessible regions (Fisher test; p=0.03; Figure 5C).

We next took an unbiased approach to identify transcription factor binding motifs that are enriched in DARs compared to all accessible regions using MEME-ChIP software (see Materials and methods). We found two motifs to be enriched in DARs compared to all regions (Figure 5D). Motif 1 (p=2×10⁻²) is recognized by HTF4 and TFE2, both of which are non-response genes in our system. Motif 2 (p=6.2×10⁻⁴²) is posited to be recognized by ZN770, E2F3, and E2F4. Both ZN770 and E2F4 are response genes in our system. DARs arise between the hypoxia and re-oxygenation conditions, and E2F4 expression increases following re-oxygenation, suggesting that it may be involved in the response. To test this hypothesis, we obtained a published ChIP-seq data set for E2F4 (Lee et al., 2011), and overlapped the 16,245 E2F4-bound regions with DARs and constitutively accessible regions. E2F4-bound regions are significantly enriched in DARs compared to constitutively accessible regions (Chi-squared test; p<2.2×10⁻¹⁶; Figure 5D). E2F4 is important for survival following ischemia in neurons, and has been suggested to be an anti-apoptotic factor in cardiomyocytes (Dingar et al., 2012; Iyirhiaro et al., 2014).

To identify additional sequence features that associate with DARs, we asked whether transposable elements (TE), a potential source of regulatory sequence subjected to chromatin-level regulation, are enriched in these sites (Du et al., 2016; He et al., 2019). We found that while three of the main TE classes, LINEs, LTRs, and DNA elements are similarly enriched in DARs compared to constitutively accessible regions; SINEs are specifically enriched in DARs (p=6.5×10⁻⁶; Figure 5E). There is an enrichment of both Alu and MIR SINE family members in DARs (p=3×10⁻⁵ and p=0.006). AluS elements, and the AluSq and AluSp sub-families, are particularly enriched within the Alu family (p=0.007 and p=0.02 respectively). A different cellular stress, heat shock, has previously been shown to remodel chromatin accessibility at Alu elements in cervical cancer cells (Kim et al., 2001).

As Alu and MIR TE sequences are notably CpG dense (Medstrand et al., 2002), we next asked about the enrichment of CpG-dense CpG islands (CGIs) in our differentially accessible regions. We found that CpG islands are enriched in DARs compared to constitutively accessible regions, whether these regions fall within 1 kb of TSS, which are typically enriched for CGIs, or not (p<2.2×10⁻¹⁶; Figure 5F).

DNA methylation state at stress-responsive genes and chromatin regions

Genes with CGI promoters are thought to allow flexibility in TSS choice compared to genes without CGI promoters (Carninci et al., 2006), and to allow for the rapid induction of gene expression in response to stimuli (Ramirez-Carrozzi et al., 2009). We therefore asked whether this promoter feature is enriched in the stress response genes. Indeed, we find that response genes are more likely to have CGI promoters than non-response genes (Chi-squared test; p=0.002).

Given the enrichment of CpG islands in gene promoters and chromatin regions that are responsive to stress, we asked whether this feature corresponded to differences in CpG DNA methylation levels in these same regions. We measured global DNA methylation levels at 766,658 CpG sites in all conditions from 13 of our individuals (no data was collected from two of the 15 individuals), together with 23 replicate samples from three individuals (see Materials and methods; Supplementary file 1). We found the expected bimodal distribution of DNA methylation Beta-values across CpGs (Beta-values represent the ratio of intensities between the methylated and unmethylated alleles; Figure 6—figure supplement 1A). Additional analyses indicated the data to be of good quality (Figure 6—figure supplements 1–2).

To determine whether steady-state DNA methylation levels mark genes or regions that will change their expression level in response to stress, we investigated baseline DNA methylation levels in the promoters of genes classified as response genes and non-response genes, as well as accessible regions and DARs. To do so, we assessed the DNA methylation level at CpGs within 200 bp upstream of the TSS in the baseline condition (condition A). The majority of the assayed CpGs were hypomethylated with a median Beta-value of less than 0.2 across genes and regions (Figure 6A). While there is no difference in median DNA methylation levels between response and non-response genes, we found that the median DNA methylation level is lower in DARs compared to all accessible regions (p<2.2×10⁻¹⁶). These data suggest that responsive chromatin regions may have specific epigenetic profiles which poise them for rapid response to stress.

Figure 6 with 4 supplements see all

Download asset Open asset

Dynamic eQTLs associate with traits and disease

Finally, we wanted to determine whether any of the dynamic eQTL SNPs or genes that we identified might also be associated with complex traits or disease. To maximize the number of potentially phenotypically-relevant genomic loci, we performed three independent analyses. We first took an unbiased, SNP-based approach by searching within a catalog of genetic variants associated with a variety of traits assayed in thousands of GWAS for overlap with our dynamic eSNPs (Buniello et al., 2019). By intersecting these two data sets, we found one induced dynamic eSNP (rs8053350) that is also associated with a measured phenotype in GWAS – varicose veins (Supplementary file 1; Figure 7—figure supplement 1A-B).

Given that our eQTL data are from the Yoruba population from West Africa and that most GWAS studies are performed in European populations, which have an inherently different LD structure, we may not expect to see strong overlap at the level of individual variants. We therefore next took an orthogonal gene-based approach, using the same GWAS catalog, to specifically investigate the genes implicated in three phenotypes associated with cardiovascular function or response to oxygen deprivation: MI, heart failure, and stroke (see Materials and methods). If the expression of these CVD genes is found to vary across individuals following hypoxia, it suggests that they may be important for mediating disease risk. Indeed, we found that six of our dynamic eGenes are implicated in these three disease states by GWAS gene mapping approaches (Supplementary file 1). This list includes the DNA damage and apoptosis factor ZC3HC1, which is implicated in MI and stroke (Figure 7—figure supplement 1C–D). Importantly, ZC3HC1 is not an eGene in LV or AA, but the SNP-gene pair is an eQTL in other tissues.

Lastly, we performed an in-depth SNP-based analysis using full summary statistics for coronary artery disease (CAD) and MI from the CARDIoGRAMplusC4D Consortium (Stitziel et al., 2016; Nikpay et al., 2015; Webb et al., 2017). We tested whether the lead eSNP from our set of dynamic eQTLs is also associated with CAD or MI and identified two loci (Supplementary file 1). This set includes the eSNP (rs12588981) for the dynamic eQTL EIF2B2 which is nominally associated with CAD (GWAS association p=4.7×10⁻⁵) and is in high and statistically significant LD with the lead GWAS variant (GWAS association p=9.9×10⁻⁸) at this locus (rs3832966; R² = 0.96, p<0.001; Figure 7; Figure 7—figure supplement 2). Given that we have incomplete power to detect dynamic eQTLs, we also used the same GWAS dataset to investigate whether any of the significant eQTL SNPs in any condition are associated with CAD or MI. We determined that eSNP (rs8105092), which regulates CARM1 in the hypoxic condition (condition B), is nominally associated with MI (GWAS association p=1.79×10⁻⁵) and in moderate, yet significant LD, with the lead GWAS variant (GWAS association p=2.8×10⁻⁹) at this locus (rs4804142; R² = 0.32, p<0.001; Figure 7—figure supplement 3), and is not an eGene in LV or AA. A handful of the eQTLs that we highlight here were tested for association with gene expression by the GTEx consortium, yet are not eQTLs in heart tissue. These results suggest that perturbation studies in relevant cell types can give insight into the molecular basis for the genetic association with complex traits and disease, which might not be gleaned from the study of post-mortem tissues.

Figure 7 with 4 supplements see all

Download asset Open asset

Studying gene expression across individuals in response to stress can reveal latent effects of genetic variation, which may contribute to higher-order phenotypes and disease. In order to understand the effects of genetic variation in a disease-relevant cell type and a disease-relevant process, we differentiated cardiomyocytes from a panel of genotyped individuals, and subjected them to hypoxia and re-oxygenation. We found hundreds of eQTLs that are revealed or suppressed following hypoxic stress (dynamic eQTLs), several of which have been associated with phenotypes measured in GWAS.

Sequencing data have been deposited in GEO under accession codes GSE144426.

1. 1. Ward MC
   2. Banovich NE
   3. Sarkar A
   4. Stephens M
   5. Gilad Y

   (2021) NCBI Gene Expression Omnibus

   ID GSE144426. Dynamic effects of genetic variation on gene expression revealed following hypoxic stress in cardiomyocytes.

   https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE144426

1. 1. GTEx Consortium

   (2020) GTEx portal V8 release

   ID V8. GTEx eQTLs.

   https://www.gtexportal.org/home/
2. 1. ENCODE Consortium

   (2012) ENCODE

   ID ENCBS684IAD. Histone ChIP-seq.

   https://www.encodeproject.org

1. 1. Afgan E
   2. Baker D
   3. Batut B
   4. van den Beek M
   5. Bouvier D
   6. Cech M
   7. Chilton J
   8. Clements D
   9. Coraor N
   10. Grüning BA
   11. Guerler A
   12. Hillman-Jackson J
   13. Hiltemann S
   14. Jalili V
   15. Rasche H
   16. Soranzo N
   17. Goecks J
   18. Taylor J
   19. Nekrutenko A
   20. Blankenberg D

   (2018) The galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update

   Nucleic Acids Research 46:W537–W544.

   https://doi.org/10.1093/nar/gky379

   • PubMed
   • Google Scholar
2. 1. Alasoo K
   2. Rodrigues J
   3. Mukhopadhyay S
   4. Knights AJ
   5. Mann AL
   6. Kundu K
   7. Hale C
   8. Dougan G
   9. Gaffney DJ
   10. HIPSCI Consortium

   (2018) Shared genetic effects on chromatin and gene expression indicate a role for enhancer priming in immune response

   Nature Genetics 50:424–431.

   https://doi.org/10.1038/s41588-018-0046-7

   • PubMed
   • Google Scholar
3. 1. Bailey TL
   2. Boden M
   3. Buske FA
   4. Frith M
   5. Grant CE
   6. Clementi L
   7. Ren J
   8. Li WW
   9. Noble WS

   (2009) MEME SUITE: tools for motif discovery and searching

   Nucleic Acids Research 37:W202–W208.

   https://doi.org/10.1093/nar/gkp335

   • Google Scholar
4. 1. Banovich NE
   2. Li YI
   3. Raj A
   4. Ward MC
   5. Greenside P
   6. Calderon D
   7. Tung PY
   8. Burnett JE
   9. Myrthil M
   10. Thomas SM
   11. Burrows CK
   12. Romero IG
   13. Pavlovic BJ
   14. Kundaje A
   15. Pritchard JK
   16. Gilad Y

   (2018) Impact of regulatory variation across human iPSCs and differentiated cells

   Genome Research 28:122–131.

   https://doi.org/10.1101/gr.224436.117

   • PubMed
   • Google Scholar
5. 1. Barreiro LB
   2. Tailleux L
   3. Pai AA
   4. Gicquel B
   5. Marioni JC
   6. Gilad Y

   (2012) Deciphering the genetic architecture of variation in the immune response to Mycobacterium tuberculosis infection

   PNAS 109:1204–1209.

   https://doi.org/10.1073/pnas.1115761109

   • PubMed
   • Google Scholar
6. 1. Batie M
   2. Frost J
   3. Frost M
   4. Wilson JW
   5. Schofield P
   6. Rocha S

   (2019) Hypoxia induces rapid changes to histone methylation and reprograms chromatin

   Science 363:1222–1226.

   https://doi.org/10.1126/science.aau5870

   • PubMed
   • Google Scholar
7. 1. Battle A
   2. Mostafavi S
   3. Zhu X
   4. Potash JB
   5. Weissman MM
   6. McCormick C
   7. Haudenschild CD
   8. Beckman KB
   9. Shi J
   10. Mei R
   11. Urban AE
   12. Montgomery SB
   13. Levinson DF
   14. Koller D

   (2014) Characterizing the genetic basis of transcriptome diversity through RNA-sequencing of 922 individuals

   Genome Research 24:14–24.

   https://doi.org/10.1101/gr.155192.113

   • PubMed
   • Google Scholar
8. 1. Benaglio P
   2. D'Antonio-Chronowska A
   3. Ma W
   4. Yang F
   5. Young Greenwald WW
   6. Donovan MKR
   7. DeBoever C
   8. Li H
   9. Drees F
   10. Singhal S
   11. Matsui H
   12. van Setten J
   13. Sotoodehnia N
   14. Gaulton KJ
   15. Smith EN
   16. D'Antonio M
   17. Rosenfeld MG
   18. Frazer KA

   (2019) Allele-specific NKX2-5 binding underlies multiple genetic associations with human electrocardiographic traits

   Nature Genetics 51:1506–1517.

   https://doi.org/10.1038/s41588-019-0499-3

   • PubMed
   • Google Scholar
9. 1. Benjamini Y
   2. Hochberg Y

   (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing

   Journal of the Royal Statistical Society: Series B 57:289–300.

   https://doi.org/10.1111/j.2517-6161.1995.tb02031.x

   • Google Scholar
10. 1. Bogdan L
    2. Barreiro L
    3. Bourque G

    (2020) Transposable elements have contributed human regulatory regions that are activated upon bacterial infection

    Philosophical Transactions of the Royal Society B: Biological Sciences 375:20190332.

    https://doi.org/10.1098/rstb.2019.0332

    • Google Scholar
11. 1. Brahimi-Horn MC
    2. Pouysségur J

    (2007) Oxygen, a source of life and stress

    FEBS Letters 581:3582–3591.

    https://doi.org/10.1016/j.febslet.2007.06.018

    • PubMed
    • Google Scholar
12. 1. Brodehl A
    2. Ebbinghaus H
    3. Deutsch MA
    4. Gummert J
    5. Gärtner A
    6. Ratnavadivel S
    7. Milting H

    (2019) Human induced pluripotent Stem-Cell-Derived cardiomyocytes as models for genetic cardiomyopathies

    International Journal of Molecular Sciences 20:4381.

    https://doi.org/10.3390/ijms20184381

    • PubMed
    • Google Scholar
13. 1. Buenrostro JD
    2. Giresi PG
    3. Zaba LC
    4. Chang HY
    5. Greenleaf WJ

    (2013) Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position

    Nature Methods 10:1213–1218.

    https://doi.org/10.1038/nmeth.2688

    • PubMed
    • Google Scholar
14. 1. Buniello A
    2. MacArthur JAL
    3. Cerezo M
    4. Harris LW
    5. Hayhurst J
    6. Malangone C
    7. McMahon A
    8. Morales J
    9. Mountjoy E
    10. Sollis E
    11. Suveges D
    12. Vrousgou O
    13. Whetzel PL
    14. Amode R
    15. Guillen JA
    16. Riat HS
    17. Trevanion SJ
    18. Hall P
    19. Junkins H
    20. Flicek P
    21. Burdett T
    22. Hindorff LA
    23. Cunningham F
    24. Parkinson H

    (2019) The NHGRI-EBI GWAS catalog of published genome-wide association studies, targeted arrays and summary statistics 2019

    Nucleic Acids Research 47:D1005–D1012.

    https://doi.org/10.1093/nar/gky1120

    • PubMed
    • Google Scholar
15. 1. Burridge PW
    2. Li YF
    3. Matsa E
    4. Wu H
    5. Ong SG
    6. Sharma A
    7. Holmström A
    8. Chang AC
    9. Coronado MJ
    10. Ebert AD
    11. Knowles JW
    12. Telli ML
    13. Witteles RM
    14. Blau HM
    15. Bernstein D
    16. Altman RB
    17. Wu JC

    (2016) Human induced pluripotent stem cell-derived cardiomyocytes recapitulate the predilection of breast Cancer patients to doxorubicin-induced cardiotoxicity

    Nature Medicine 22:547–556.

    https://doi.org/10.1038/nm.4087

    • PubMed
    • Google Scholar
16. 1. Calderon D
    2. Nguyen MLT
    3. Mezger A
    4. Kathiria A
    5. Müller F
    6. Nguyen V
    7. Lescano N
    8. Wu B
    9. Trombetta J
    10. Ribado JV
    11. Knowles DA
    12. Gao Z
    13. Blaeschke F
    14. Parent AV
    15. Burt TD
    16. Anderson MS
    17. Criswell LA
    18. Greenleaf WJ
    19. Marson A
    20. Pritchard JK

    (2019) Landscape of stimulation-responsive chromatin across diverse human immune cells

    Nature Genetics 51:1494–1505.

    https://doi.org/10.1038/s41588-019-0505-9

    • PubMed
    • Google Scholar
17. 1. Çalışkan M
    2. Baker SW
    3. Gilad Y
    4. Ober C

    (2015) Host genetic variation influences gene expression response to rhinovirus infection

    PLOS Genetics 11:e1005111.

    https://doi.org/10.1371/journal.pgen.1005111

    • PubMed
    • Google Scholar
18. 1. Carninci P
    2. Sandelin A
    3. Lenhard B
    4. Katayama S
    5. Shimokawa K
    6. Ponjavic J
    7. Semple CA
    8. Taylor MS
    9. Engström PG
    10. Frith MC
    11. Forrest AR
    12. Alkema WB
    13. Tan SL
    14. Plessy C
    15. Kodzius R
    16. Ravasi T
    17. Kasukawa T
    18. Fukuda S
    19. Kanamori-Katayama M
    20. Kitazume Y
    21. Kawaji H
    22. Kai C
    23. Nakamura M
    24. Konno H
    25. Nakano K
    26. Mottagui-Tabar S
    27. Arner P
    28. Chesi A
    29. Gustincich S
    30. Persichetti F
    31. Suzuki H
    32. Grimmond SM
    33. Wells CA
    34. Orlando V
    35. Wahlestedt C
    36. Liu ET
    37. Harbers M
    38. Kawai J
    39. Bajic VB
    40. Hume DA
    41. Hayashizaki Y

    (2006) Genome-wide analysis of mammalian promoter architecture and evolution

    Nature Genetics 38:626–635.

    https://doi.org/10.1038/ng1789

    • PubMed
    • Google Scholar
19. 1. Carreau A
    2. El Hafny-Rahbi B
    3. Matejuk A
    4. Grillon C
    5. Kieda C

    (2011) Why is the partial oxygen pressure of human tissues a crucial parameter? small molecules and hypoxia

    Journal of Cellular and Molecular Medicine 15:1239–1253.

    https://doi.org/10.1111/j.1582-4934.2011.01258.x

    • PubMed
    • Google Scholar
20. 1. Chang SL
    2. Huang YL
    3. Lee MC
    4. Hu S
    5. Hsiao YC
    6. Chang SW
    7. Chang CJ
    8. Chen PC

    (2018) Association of varicose veins with incident venous thromboembolism and peripheral artery disease

    Jama 319:807–817.

    https://doi.org/10.1001/jama.2018.0246

    • PubMed
    • Google Scholar
21. 1. Cuomo ASE
    2. Seaton DD
    3. McCarthy DJ
    4. Martinez I
    5. Bonder MJ
    6. Garcia-Bernardo J
    7. Amatya S
    8. Madrigal P
    9. Isaacson A
    10. Buettner F
    11. Knights A
    12. Natarajan KN
    13. Vallier L
    14. Marioni JC
    15. Chhatriwala M
    16. Stegle O
    17. HipSci Consortium

    (2020) Single-cell RNA-sequencing of differentiating iPS cells reveals dynamic genetic effects on gene expression

    Nature Communications 11:810.

    https://doi.org/10.1038/s41467-020-14457-z

    • PubMed
    • Google Scholar
22. 1. Davis CA
    2. Hitz BC
    3. Sloan CA
    4. Chan ET
    5. Davidson JM
    6. Gabdank I
    7. Hilton JA
    8. Jain K
    9. Baymuradov UK
    10. Narayanan AK
    11. Onate KC
    12. Graham K
    13. Miyasato SR
    14. Dreszer TR
    15. Strattan JS
    16. Jolanki O
    17. Tanaka FY
    18. Cherry JM

    (2018) The encyclopedia of DNA elements (ENCODE): data portal update

    Nucleic Acids Research 46:D794–D801.

    https://doi.org/10.1093/nar/gkx1081

    • PubMed
    • Google Scholar
23. 1. de la Roche J
    2. Angsutararux P
    3. Kempf H
    4. Janan M
    5. Bolesani E
    6. Thiemann S
    7. Wojciechowski D
    8. Coffee M
    9. Franke A
    10. Schwanke K
    11. Leffler A
    12. Luanpitpong S
    13. Issaragrisil S
    14. Fischer M
    15. Zweigerdt R

    (2019) Comparing human iPSC-cardiomyocytes versus HEK293T cells unveils disease-causing effects of brugada mutation A735V of NaV1.5 sodium channels

    Scientific Reports 9:11173.

    https://doi.org/10.1038/s41598-019-47632-4

    • Google Scholar
24. 1. Delaneau O
    2. Ongen H
    3. Brown AA
    4. Fort A
    5. Panousis NI
    6. Dermitzakis ET

    (2017) A complete tool set for molecular QTL discovery and analysis

    Nature Communications 8:15452.

    https://doi.org/10.1038/ncomms15452

    • PubMed
    • Google Scholar
25. 1. Dingar D
    2. Konecny F
    3. Zou J
    4. Sun X
    5. von Harsdorf R

    (2012) Anti-apoptotic function of the E2F transcription factor 4 (E2F4)/p130, a member of retinoblastoma gene family in cardiac myocytes

    Journal of Molecular and Cellular Cardiology 53:820–828.

    https://doi.org/10.1016/j.yjmcc.2012.09.004

    • Google Scholar
26. 1. Dombroski BA
    2. Nayak RR
    3. Ewens KG
    4. Ankener W
    5. Cheung VG
    6. Spielman RS

    (2010) Gene expression and genetic variation in response to endoplasmic reticulum stress in human cells

    The American Journal of Human Genetics 86:719–729.

    https://doi.org/10.1016/j.ajhg.2010.03.017

    • PubMed
    • Google Scholar
27. Preprint

    1. Donovan M-C
    2. D'Antonio A.;
    3. M.; Frazer KA

    (2019) Cellular deconvolution of GTEx tissues powers eQTL studies to discover thousands of novel disease and cell-type associated regulatory variants

    bioRxiv.

    https://doi.org/10.1101/671040

    • Google Scholar
28. 1. Du P
    2. Kibbe WA
    3. Lin SM

    (2008) Lumi: a pipeline for processing illumina microarray

    Bioinformatics 24:1547–1548.

    https://doi.org/10.1093/bioinformatics/btn224

    • PubMed
    • Google Scholar
29. 1. Du J
    2. Leung A
    3. Trac C
    4. Lee M
    5. Parks BW
    6. Lusis AJ
    7. Natarajan R
    8. Schones DE

    (2016) Chromatin variation associated with liver metabolism is mediated by transposable elements

    Epigenetics & Chromatin 9:28.

    https://doi.org/10.1186/s13072-016-0078-0

    • Google Scholar
30. 1. Dzau VJ
    2. Antman EM
    3. Black HR
    4. Hayes DL
    5. Manson JE
    6. Plutzky J
    7. Popma JJ
    8. Stevenson W

    (2006) The cardiovascular disease continuum validated: clinical evidence of improved patient outcomes: part I: pathophysiology and clinical trial evidence (risk factors through stable coronary artery disease)

    Circulation 114:2850–2870.

    https://doi.org/10.1161/CIRCULATIONAHA.106.655688

    • PubMed
    • Google Scholar
31. 1. ENCODE Project Consortium

    (2012) An integrated encyclopedia of DNA elements in the human genome

    Nature 489:57–74.

    https://doi.org/10.1038/nature11247

    • PubMed
    • Google Scholar
32. 1. Feige E
    2. Yokoyama S
    3. Levy C
    4. Khaled M
    5. Igras V
    6. Lin RJ
    7. Lee S
    8. Widlund HR
    9. Granter SR
    10. Kung AL
    11. Fisher DE

    (2011) Hypoxia-induced transcriptional repression of the melanoma-associated oncogene MITF

    PNAS 108:E924–E933.

    https://doi.org/10.1073/pnas.1106351108

    • PubMed
    • Google Scholar
33. 1. Ferrari R
    2. de Llobet Cucalon LI
    3. Di Vona C
    4. Le Dilly F
    5. Vidal E
    6. Lioutas A
    7. Oliete JQ
    8. Jochem L
    9. Cutts E
    10. Dieci G
    11. Vannini A
    12. Teichmann M
    13. de la Luna S
    14. Beato M

    (2020) TFIIIC binding to alu elements controls gene expression via chromatin looping and histone acetylation

    Molecular Cell 77:475–487.

    https://doi.org/10.1016/j.molcel.2019.10.020

    • PubMed
    • Google Scholar
34. 1. Findley AS
    2. Richards AL
    3. Petrini C
    4. Alazizi A
    5. Doman E
    6. Shanku AG
    7. Davis GO
    8. Hauff N
    9. Sorokin Y
    10. Wen X
    11. Pique-Regi R
    12. Luca F

    (2019) Interpreting coronary artery disease risk through Gene-Environment interactions in gene regulation

    Genetics 213:651–663.

    https://doi.org/10.1534/genetics.119.302419

    • PubMed
    • Google Scholar
35. 1. Fukaya E
    2. Flores AM
    3. Lindholm D
    4. Gustafsson S
    5. Zanetti D
    6. Ingelsson E
    7. Leeper NJ

    (2018) Clinical and genetic determinants of varicose veins

    Circulation 138:2869–2880.

    https://doi.org/10.1161/CIRCULATIONAHA.118.035584

    • PubMed
    • Google Scholar
36. 1. Ghorbel MT
    2. Cherif M
    3. Jenkins E
    4. Mokhtari A
    5. Kenny D
    6. Angelini GD
    7. Caputo M

    (2010) Transcriptomic analysis of patients with tetralogy of fallot reveals the effect of chronic hypoxia on myocardial gene expression

    The Journal of Thoracic and Cardiovascular Surgery 140:337–345.

    https://doi.org/10.1016/j.jtcvs.2009.12.055

    • PubMed
    • Google Scholar
37. 1. Giordano FJ

    (2005) Oxygen, oxidative stress, hypoxia, and heart failure

    Journal of Clinical Investigation 115:500–508.

    https://doi.org/10.1172/JCI200524408

    • PubMed
    • Google Scholar
38. 1. GTEx Consortium
    2. Laboratory, Data Analysis &Coordinating Center (LDACC)—Analysis Working Group
    3. Statistical Methods groups—Analysis Working Group
    4. Enhancing GTEx (eGTEx) groups
    5. NIH Common Fund
    6. NIH/NCI
    7. NIH/NHGRI
    8. NIH/NIMH
    9. NIH/NIDA
    10. Biospecimen Collection Source Site—NDRI
    11. Biospecimen Collection Source Site—RPCI
    12. Biospecimen Core Resource—VARI
    13. Brain Bank Repository—University of Miami Brain Endowment Bank
    14. Leidos Biomedical—Project Management
    15. ELSI Study
    16. Genome Browser Data Integration &Visualization—EBI
    17. Genome Browser Data Integration &Visualization—UCSC Genomics Institute, University of California Santa Cruz
    18. Lead analysts:
    19. Laboratory, Data Analysis &Coordinating Center (LDACC):
    20. NIH program management:
    21. Biospecimen collection:
    22. Pathology:
    23. eQTL manuscript working group:
    24. Battle A
    25. Brown CD
    26. Engelhardt BE
    27. Montgomery SB

    (2017) Genetic effects on gene expression across human tissues

    Nature 550:204–213.

    https://doi.org/10.1038/nature24277

    • PubMed
    • Google Scholar
39. 1. Gusev A
    2. Ko A
    3. Shi H
    4. Bhatia G
    5. Chung W
    6. Penninx BW
    7. Jansen R
    8. de Geus EJ
    9. Boomsma DI
    10. Wright FA
    11. Sullivan PF
    12. Nikkola E
    13. Alvarez M
    14. Civelek M
    15. Lusis AJ
    16. Lehtimäki T
    17. Raitoharju E
    18. Kähönen M
    19. Seppälä I
    20. Raitakari OT
    21. Kuusisto J
    22. Laakso M
    23. Price AL
    24. Pajukanta P
    25. Pasaniuc B

    (2016) Integrative approaches for large-scale transcriptome-wide association studies

    Nature Genetics 48:245–252.

    https://doi.org/10.1038/ng.3506

    • PubMed
    • Google Scholar
40. 1. Hardin M
    2. Cho MH
    3. McDonald ML
    4. Wan E
    5. Lomas DA
    6. Coxson HO
    7. MacNee W
    8. Vestbo J
    9. Yates JC
    10. Agusti A
    11. Calverley PM
    12. Celli B
    13. Crim C
    14. Rennard S
    15. Wouters E
    16. Bakke P
    17. Bhatt SP
    18. Kim V
    19. Ramsdell J
    20. Regan EA
    21. Make BJ
    22. Hokanson JE
    23. Crapo JD
    24. Beaty TH
    25. Hersh CP
    26. ECLIPSE and COPDGene Investigators
    27. COPDGene Investigators—clinical centers

    (2016) A genome-wide analysis of the response to inhaled β2-agonists in chronic obstructive pulmonary disease

    The Pharmacogenomics Journal 16:326–335.

    https://doi.org/10.1038/tpj.2015.65

    • PubMed
    • Google Scholar
41. 1. Hartley I
    2. Elkhoury FF
    3. Heon Shin J
    4. Xie B
    5. Gu X
    6. Gao Y
    7. Zhou D
    8. Haddad GG

    (2013) Long-lasting changes in DNA methylation following short-term hypoxic exposure in primary hippocampal neuronal cultures

    PLOS ONE 8:e77859.

    https://doi.org/10.1371/journal.pone.0077859

    • PubMed
    • Google Scholar
42. 1. He J
    2. Fu X
    3. Zhang M
    4. He F
    5. Li W
    6. Abdul MM
    7. Zhou J
    8. Sun L
    9. Chang C
    10. Li Y
    11. Liu H
    12. Wu K
    13. Babarinde IA
    14. Zhuang Q
    15. Loh YH
    16. Chen J
    17. Esteban MA
    18. Hutchins AP

    (2019) Transposable elements are regulated by context-specific patterns of chromatin marks in mouse embryonic stem cells

    Nature Communications 10:34.

    https://doi.org/10.1038/s41467-018-08006-y

    • PubMed
    • Google Scholar
43. 1. Heinig M
    2. Adriaens ME
    3. Schafer S
    4. van Deutekom HWM
    5. Lodder EM
    6. Ware JS
    7. Schneider V
    8. Felkin LE
    9. Creemers EE
    10. Meder B
    11. Katus HA
    12. Rühle F
    13. Stoll M
    14. Cambien F
    15. Villard E
    16. Charron P
    17. Varro A
    18. Bishopric NH
    19. George AL
    20. Dos Remedios C
    21. Moreno-Moral A
    22. Pesce F
    23. Bauerfeind A
    24. Rüschendorf F
    25. Rintisch C
    26. Petretto E
    27. Barton PJ
    28. Cook SA
    29. Pinto YM
    30. Bezzina CR
    31. Hubner N

    (2017) Natural genetic variation of the cardiac transcriptome in non-diseased donors and patients with dilated cardiomyopathy

    Genome Biology 18:170.

    https://doi.org/10.1186/s13059-017-1286-z

    • PubMed
    • Google Scholar
44. 1. Huang daW
    2. Sherman BT
    3. Lempicki RA

    (2009a) Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists

    Nucleic Acids Research 37:1–13.

    https://doi.org/10.1093/nar/gkn923

    • PubMed
    • Google Scholar
45. 1. Huang daW
    2. Sherman BT
    3. Lempicki RA

    (2009b) Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources

    Nature Protocols 4:44–57.

    https://doi.org/10.1038/nprot.2008.211

    • PubMed
    • Google Scholar
46. 1. IBC 50K CAD Consortium

    (2011) Large-scale gene-centric analysis identifies novel variants for coronary artery disease

    PLOS Genetics 7:e1002260.

    https://doi.org/10.1371/journal.pgen.1002260

    • PubMed
    • Google Scholar
47. 1. Iyirhiaro GO
    2. Zhang Y
    3. Estey C
    4. O'Hare MJ
    5. Safarpour F
    6. Parsanejad M
    7. Wang S
    8. Abdel-Messih E
    9. Callaghan SM
    10. During MJ
    11. Slack RS
    12. Park DS

    (2014) Regulation of ischemic neuronal death by E2F4-p130 protein complexes

    Journal of Biological Chemistry 289:18202–18213.

    https://doi.org/10.1074/jbc.M114.574145

    • PubMed
    • Google Scholar
48. 1. Jagannathan L
    2. Cuddapah S
    3. Costa M

    (2016) Oxidative stress under ambient and physiological oxygen tension in tissue culture

    Current Pharmacology Reports 2:64–72.

    https://doi.org/10.1007/s40495-016-0050-5

    • PubMed
    • Google Scholar
49. 1. Joehanes R
    2. Zhang X
    3. Huan T
    4. Yao C
    5. Ying SX
    6. Nguyen QT
    7. Demirkale CY
    8. Feolo ML
    9. Sharopova NR
    10. Sturcke A
    11. Schäffer AA
    12. Heard-Costa N
    13. Chen H
    14. Liu PC
    15. Wang R
    16. Woodhouse KA
    17. Tanriverdi K
    18. Freedman JE
    19. Raghavachari N
    20. Dupuis J
    21. Johnson AD
    22. O'Donnell CJ
    23. Levy D
    24. Munson PJ

    (2017) Integrated genome-wide analysis of expression quantitative trait loci aids interpretation of genomic association studies

    Genome Biology 18:16.

    https://doi.org/10.1186/s13059-016-1142-6

    • PubMed
    • Google Scholar
50. 1. Jurka J

    (2000) Repbase update: a database and an electronic

    Journal of Repetitive Elements Trends Genet 16:418–420.

    https://doi.org/10.1016/S0168-9525(00)02093-X

    • Google Scholar
51. 1. Kariuki SN
    2. Maranville JC
    3. Baxter SS
    4. Jeong C
    5. Nakagome S
    6. Hrusch CL
    7. Witonsky DB
    8. Sperling AI
    9. Di Rienzo A

    (2016) Mapping variation in cellular and transcriptional response to 1,25-Dihydroxyvitamin D3 in peripheral blood mononuclear cells

    PLOS ONE 11:e0159779.

    https://doi.org/10.1371/journal.pone.0159779

    • PubMed
    • Google Scholar
52. 1. Karolchik D
    2. Hinrichs AS
    3. Furey TS
    4. Roskin KM
    5. Sugnet CW
    6. Haussler D
    7. Kent WJ

    (2004) The UCSC table browser data retrieval tool

    Nucleic Acids Research 32:493–496.

    https://doi.org/10.1093/nar/gkh103

    • PubMed
    • Google Scholar
53. 1. Kim C
    2. Rubin CM
    3. Schmid CW

    (2001) Genome-wide chromatin remodeling modulates the alu heat shock response

    Gene 276:127–133.

    https://doi.org/10.1016/S0378-1119(01)00639-4

    • PubMed
    • Google Scholar
54. 1. Kim-Hellmuth S
    2. Bechheim M
    3. Pütz B
    4. Mohammadi P
    5. Nédélec Y
    6. Giangreco N
    7. Becker J
    8. Kaiser V
    9. Fricker N
    10. Beier E
    11. Boor P
    12. Castel SE
    13. Nöthen MM
    14. Barreiro LB
    15. Pickrell JK
    16. Müller-Myhsok B
    17. Lappalainen T
    18. Schumacher J
    19. Hornung V

    (2017) Genetic regulatory effects modified by immune activation contribute to autoimmune disease associations

    Nature Communications 8:266.

    https://doi.org/10.1038/s41467-017-00366-1

    • PubMed
    • Google Scholar
55. 1. Knowles DA
    2. Burrows CK
    3. Blischak JD
    4. Patterson KM
    5. Serie DJ
    6. Norton N
    7. Ober C
    8. Pritchard JK
    9. Gilad Y

    (2018) Determining the genetic basis of anthracycline-cardiotoxicity by molecular response QTL mapping in induced cardiomyocytes

    eLife 7:e33480.

    https://doi.org/10.7554/eLife.33480

    • PubMed
    • Google Scholar
56. 1. Koopmann TT
    2. Adriaens ME
    3. Moerland PD
    4. Marsman RF
    5. Westerveld ML
    6. Lal S
    7. Zhang T
    8. Simmons CQ
    9. Baczko I
    10. dos Remedios C
    11. Bishopric NH
    12. Varro A
    13. George AL
    14. Lodder EM
    15. Bezzina CR

    (2014) Genome-wide identification of expression quantitative trait loci (eQTLs) in human heart

    PLOS ONE 9:e97380.

    https://doi.org/10.1371/journal.pone.0097380

    • PubMed
    • Google Scholar
57. 1. Lambert SA
    2. Jolma A
    3. Campitelli LF
    4. Das PK
    5. Yin Y
    6. Albu M
    7. Chen X
    8. Taipale J
    9. Hughes TR
    10. Weirauch MT

    (2018) The human transcription factors

    Cell 175:598–599.

    https://doi.org/10.1016/j.cell.2018.09.045

    • Google Scholar
58. 1. Langmead B
    2. Salzberg SL

    (2012) Fast gapped-read alignment with bowtie 2

    Nature Methods 9:357–359.

    https://doi.org/10.1038/nmeth.1923

    • PubMed
    • Google Scholar
59. 1. Law CW
    2. Chen Y
    3. Shi W
    4. Smyth GK

    (2014) Voom: precision weights unlock linear model analysis tools for RNA-seq read counts

    Genome Biology 15:R29.

    https://doi.org/10.1186/gb-2014-15-2-r29

    • PubMed
    • Google Scholar
60. 1. Lee BK
    2. Bhinge AA
    3. Iyer VR

    (2011) Wide-ranging functions of E2F4 in transcriptional activation and repression revealed by genome-wide analysis

    Nucleic Acids Research 39:3558–3573.

    https://doi.org/10.1093/nar/gkq1313

    • PubMed
    • Google Scholar
61. 1. Liao Y
    2. Smyth GK
    3. Shi W

    (2013) The subread aligner: fast, accurate and scalable read mapping by seed-and-vote

    Nucleic Acids Research 41:e108.

    https://doi.org/10.1093/nar/gkt214

    • PubMed
    • Google Scholar
62. 1. Liao Y
    2. Smyth GK
    3. Shi W

    (2014) featureCounts: an efficient general purpose program for assigning sequence reads to genomic features

    Bioinformatics 30:923–930.

    https://doi.org/10.1093/bioinformatics/btt656

    • PubMed
    • Google Scholar
63. 1. Lin H
    2. Dolmatova EV
    3. Morley MP
    4. Lunetta KL
    5. McManus DD
    6. Magnani JW
    7. Margulies KB
    8. Hakonarson H
    9. del Monte F
    10. Benjamin EJ
    11. Cappola TP
    12. Ellinor PT

    (2014) Gene expression and genetic variation in human atria

    Heart Rhythm 11:266–271.

    https://doi.org/10.1016/j.hrthm.2013.10.051

    • PubMed
    • Google Scholar
64. 1. Ma N
    2. Zhang JZ
    3. Itzhaki I
    4. Zhang SL
    5. Chen H
    6. Haddad F
    7. Kitani T
    8. Wilson KD
    9. Tian L
    10. Shrestha R
    11. Wu H
    12. Lam CK
    13. Sayed N
    14. Wu JC

    (2018) Determining the pathogenicity of a genomic variant of uncertain significance using CRISPR/Cas9 and Human-Induced pluripotent stem cells

    Circulation 138:2666–2681.

    https://doi.org/10.1161/CIRCULATIONAHA.117.032273

    • PubMed
    • Google Scholar
65. 1. Machanick P
    2. Bailey TL

    (2011) MEME-ChIP: motif analysis of large DNA datasets

    Bioinformatics 27:1696–1697.

    https://doi.org/10.1093/bioinformatics/btr189

    • PubMed
    • Google Scholar
66. 1. Machiela MJ
    2. Chanock SJ

    (2015) LDlink: a web-based application for exploring population-specific haplotype structure and linking correlated alleles of possible functional variants

    Bioinformatics 31:3555–3557.

    https://doi.org/10.1093/bioinformatics/btv402

    • PubMed
    • Google Scholar
67. 1. Maeda K
    2. Tsutamoto T
    3. Wada A
    4. Hisanaga T
    5. Kinoshita M

    (1998) Plasma brain natriuretic peptide as a biochemical marker of high left ventricular end-diastolic pressure in patients with symptomatic left ventricular dysfunction

    American Heart Journal 135:825–832.

    https://doi.org/10.1016/S0002-8703(98)70041-9

    • PubMed
    • Google Scholar
68. 1. Manry J
    2. Nédélec Y
    3. Fava VM
    4. Cobat A
    5. Orlova M
    6. Thuc NV
    7. Thai VH
    8. Laval G
    9. Barreiro LB
    10. Schurr E

    (2017) Deciphering the genetic control of gene expression following Mycobacterium leprae antigen stimulation

    PLOS Genetics 13:e1006952.

    https://doi.org/10.1371/journal.pgen.1006952

    • PubMed
    • Google Scholar
69. 1. Maranville JC
    2. Luca F
    3. Richards AL
    4. Wen X
    5. Witonsky DB
    6. Baxter S
    7. Stephens M
    8. Di Rienzo A

    (2011) Interactions between glucocorticoid treatment and cis-regulatory polymorphisms contribute to cellular response phenotypes

    PLOS Genetics 7:e1002162.

    https://doi.org/10.1371/journal.pgen.1002162

    • PubMed
    • Google Scholar
70. 1. McDermott-Roe C
    2. Lv W
    3. Maximova T
    4. Wada S
    5. Bukowy J
    6. Marquez M
    7. Lai S
    8. Shehu A
    9. Benjamin I
    10. Geurts A
    11. Musunuru K

    (2019) Investigation of a dilated cardiomyopathy-associated variant in BAG3 using genome-edited iPSC-derived cardiomyocytes

    JCI Insight 4:128799.

    https://doi.org/10.1172/jci.insight.128799

    • PubMed
    • Google Scholar
71. 1. Medstrand P
    2. van de Lagemaat LN
    3. Mager DL

    (2002) Retroelement distributions in the human genome: variations associated with age and proximity to genes

    Genome Research 12:1483–1495.

    https://doi.org/10.1101/gr.388902

    • PubMed
    • Google Scholar
72. 1. Narravula S
    2. Colgan SP

    (2001) Hypoxia-inducible factor 1-mediated inhibition of peroxisome proliferator-activated receptor alpha expression during hypoxia

    The Journal of Immunology 166:7543–7548.

    https://doi.org/10.4049/jimmunol.166.12.7543

    • PubMed
    • Google Scholar
73. 1. Nédélec Y
    2. Sanz J
    3. Baharian G
    4. Szpiech ZA
    5. Pacis A
    6. Dumaine A
    7. Grenier JC
    8. Freiman A
    9. Sams AJ
    10. Hebert S
    11. Pagé Sabourin A
    12. Luca F
    13. Blekhman R
    14. Hernandez RD
    15. Pique-Regi R
    16. Tung J
    17. Yotova V
    18. Barreiro LB

    (2016) Genetic ancestry and natural selection drive population differences in immune responses to pathogens

    Cell 167:657–669.

    https://doi.org/10.1016/j.cell.2016.09.025

    • PubMed
    • Google Scholar
74. 1. Nikpay M
    2. Goel A
    3. Won HH
    4. Hall LM
    5. Willenborg C
    6. Kanoni S
    7. Saleheen D
    8. Kyriakou T
    9. Nelson CP
    10. Hopewell JC
    11. Webb TR
    12. Zeng L
    13. Dehghan A
    14. Alver M
    15. Armasu SM
    16. Auro K
    17. Bjonnes A
    18. Chasman DI
    19. Chen S
    20. Ford I
    21. Franceschini N
    22. Gieger C
    23. Grace C
    24. Gustafsson S
    25. Huang J
    26. Hwang SJ
    27. Kim YK
    28. Kleber ME
    29. Lau KW
    30. Lu X
    31. Lu Y
    32. Lyytikäinen LP
    33. Mihailov E
    34. Morrison AC
    35. Pervjakova N
    36. Qu L
    37. Rose LM
    38. Salfati E
    39. Saxena R
    40. Scholz M
    41. Smith AV
    42. Tikkanen E
    43. Uitterlinden A
    44. Yang X
    45. Zhang W
    46. Zhao W
    47. de Andrade M
    48. de Vries PS
    49. van Zuydam NR
    50. Anand SS
    51. Bertram L
    52. Beutner F
    53. Dedoussis G
    54. Frossard P
    55. Gauguier D
    56. Goodall AH
    57. Gottesman O
    58. Haber M
    59. Han BG
    60. Huang J
    61. Jalilzadeh S
    62. Kessler T
    63. König IR
    64. Lannfelt L
    65. Lieb W
    66. Lind L
    67. Lindgren CM
    68. Lokki ML
    69. Magnusson PK
    70. Mallick NH
    71. Mehra N
    72. Meitinger T
    73. Memon FU
    74. Morris AP
    75. Nieminen MS
    76. Pedersen NL
    77. Peters A
    78. Rallidis LS
    79. Rasheed A
    80. Samuel M
    81. Shah SH
    82. Sinisalo J
    83. Stirrups KE
    84. Trompet S
    85. Wang L
    86. Zaman KS
    87. Ardissino D
    88. Boerwinkle E
    89. Borecki IB
    90. Bottinger EP
    91. Buring JE
    92. Chambers JC
    93. Collins R
    94. Cupples LA
    95. Danesh J
    96. Demuth I
    97. Elosua R
    98. Epstein SE
    99. Esko T
    100. Feitosa MF
    101. Franco OH
    102. Franzosi MG
    103. Granger CB
    104. Gu D
    105. Gudnason V
    106. Hall AS
    107. Hamsten A
    108. Harris TB
    109. Hazen SL
    110. Hengstenberg C
    111. Hofman A
    112. Ingelsson E
    113. Iribarren C
    114. Jukema JW
    115. Karhunen PJ
    116. Kim BJ
    117. Kooner JS
    118. Kullo IJ
    119. Lehtimäki T
    120. Loos RJF
    121. Melander O
    122. Metspalu A
    123. März W
    124. Palmer CN
    125. Perola M
    126. Quertermous T
    127. Rader DJ
    128. Ridker PM
    129. Ripatti S
    130. Roberts R
    131. Salomaa V
    132. Sanghera DK
    133. Schwartz SM
    134. Seedorf U
    135. Stewart AF
    136. Stott DJ
    137. Thiery J
    138. Zalloua PA
    139. O'Donnell CJ
    140. Reilly MP
    141. Assimes TL
    142. Thompson JR
    143. Erdmann J
    144. Clarke R
    145. Watkins H
    146. Kathiresan S
    147. McPherson R
    148. Deloukas P
    149. Schunkert H
    150. Samani NJ
    151. Farrall M

    (2015) A comprehensive 1,000 Genomes-based genome-wide association meta-analysis of coronary artery disease

    Nature Genetics 47:1121–1130.

    https://doi.org/10.1038/ng.3396

    • PubMed
    • Google Scholar
75. 1. Ongen H
    2. Buil A
    3. Brown AA
    4. Dermitzakis ET
    5. Delaneau O

    (2016) Fast and efficient QTL mapper for thousands of molecular phenotypes

    Bioinformatics 32:1479–1485.

    https://doi.org/10.1093/bioinformatics/btv722

    • PubMed
    • Google Scholar
76. 1. Pacis A
    2. Tailleux L
    3. Morin AM
    4. Lambourne J
    5. MacIsaac JL
    6. Yotova V
    7. Dumaine A
    8. Danckaert A
    9. Luca F
    10. Grenier JC
    11. Hansen KD
    12. Gicquel B
    13. Yu M
    14. Pai A
    15. He C
    16. Tung J
    17. Pastinen T
    18. Kobor MS
    19. Pique-Regi R
    20. Gilad Y
    21. Barreiro LB

    (2015) Bacterial infection remodels the DNA methylation landscape of human dendritic cells

    Genome Research 25:1801–1811.

    https://doi.org/10.1101/gr.192005.115

    • PubMed
    • Google Scholar
77. 1. Pacis A
    2. Mailhot-Léonard F
    3. Tailleux L
    4. Randolph HE
    5. Yotova V
    6. Dumaine A
    7. Grenier JC
    8. Barreiro LB

    (2019) Gene activation precedes DNA demethylation in response to infection in human dendritic cells

    PNAS 116:6938–6943.

    https://doi.org/10.1073/pnas.1814700116

    • PubMed
    • Google Scholar
78. 1. Panopoulos AD
    2. D'Antonio M
    3. Benaglio P
    4. Williams R
    5. Hashem SI
    6. Schuldt BM
    7. DeBoever C
    8. Arias AD
    9. Garcia M
    10. Nelson BC
    11. Harismendy O
    12. Jakubosky DA
    13. Donovan MKR
    14. Greenwald WW
    15. Farnam K
    16. Cook M
    17. Borja V
    18. Miller CA
    19. Grinstein JD
    20. Drees F
    21. Okubo J
    22. Diffenderfer KE
    23. Hishida Y
    24. Modesto V
    25. Dargitz CT
    26. Feiring R
    27. Zhao C
    28. Aguirre A
    29. McGarry TJ
    30. Matsui H
    31. Li H
    32. Reyna J
    33. Rao F
    34. O'Connor DT
    35. Yeo GW
    36. Evans SM
    37. Chi NC
    38. Jepsen K
    39. Nariai N
    40. Müller FJ
    41. Goldstein LSB
    42. Izpisua Belmonte JC
    43. Adler E
    44. Loring JF
    45. Berggren WT
    46. D'Antonio-Chronowska A
    47. Smith EN
    48. Frazer KA

    (2017) iPSCORE: a resource of 222 iPSC lines enabling functional characterization of genetic variation across a variety of cell types

    Stem Cell Reports 8:1086–1100.

    https://doi.org/10.1016/j.stemcr.2017.03.012

    • PubMed
    • Google Scholar
79. 1. Pavlides JM
    2. Zhu Z
    3. Gratten J
    4. McRae AF
    5. Wray NR
    6. Yang J

    (2016) Predicting gene targets from integrative analyses of summary data from GWAS and eQTL studies for 28 human complex traits

    Genome Medicine 8:84.

    https://doi.org/10.1186/s13073-016-0338-4

    • PubMed
    • Google Scholar
80. 1. Pavlovic BJ
    2. Blake LE
    3. Roux J
    4. Chavarria C
    5. Gilad Y

    (2018) A comparative assessment of human and chimpanzee iPSC-derived cardiomyocytes with primary heart tissues

    Scientific Reports 8:15312.

    https://doi.org/10.1038/s41598-018-33478-9

    • PubMed
    • Google Scholar
81. 1. Pickrell JK
    2. Marioni JC
    3. Pai AA
    4. Degner JF
    5. Engelhardt BE
    6. Nkadori E
    7. Veyrieras J-B
    8. Stephens M
    9. Gilad Y
    10. Pritchard JK

    (2010) Understanding mechanisms underlying human gene expression variation with RNA sequencing

    Nature 464:768–772.

    https://doi.org/10.1038/nature08872

    • Google Scholar
82. 1. Pinto AR
    2. Ilinykh A
    3. Ivey MJ
    4. Kuwabara JT
    5. D'Antoni ML
    6. Debuque R
    7. Chandran A
    8. Wang L
    9. Arora K
    10. Rosenthal NA
    11. Tallquist MD

    (2016) Revisiting cardiac cellular composition

    Circulation Research 118:400–409.

    https://doi.org/10.1161/CIRCRESAHA.115.307778

    • PubMed
    • Google Scholar
83. 1. Platt JL
    2. Salama R
    3. Smythies J
    4. Choudhry H
    5. Davies JO
    6. Hughes JR
    7. Ratcliffe PJ
    8. Mole DR

    (2016) Capture-C reveals preformed chromatin interactions between HIF-binding sites and distant promoters

    EMBO Reports 17:1410–1421.

    https://doi.org/10.15252/embr.201642198

    • PubMed
    • Google Scholar
84. 1. Poetsch AR
    2. Boulton SJ
    3. Luscombe NM

    (2018) Genomic landscape of oxidative DNA damage and repair reveals regioselective protection from mutagenesis

    Genome Biology 19:215.

    https://doi.org/10.1186/s13059-018-1582-2

    • PubMed
    • Google Scholar
85. 1. Quinlan AR
    2. Hall IM

    (2010) BEDTools: a flexible suite of utilities for comparing genomic features

    Bioinformatics 26:841–842.

    https://doi.org/10.1093/bioinformatics/btq033

    • PubMed
    • Google Scholar
86. 1. Ramirez-Carrozzi VR
    2. Braas D
    3. Bhatt DM
    4. Cheng CS
    5. Hong C
    6. Doty KR
    7. Black JC
    8. Hoffmann A
    9. Carey M
    10. Smale ST

    (2009) A unifying model for the selective regulation of inducible transcription by CpG islands and nucleosome remodeling

    Cell 138:114–128.

    https://doi.org/10.1016/j.cell.2009.04.020

    • PubMed
    • Google Scholar
87. 1. Risso D
    2. Ngai J
    3. Speed TP
    4. Dudoit S

    (2014) Normalization of RNA-seq data using factor analysis of control genes or samples

    Nature Biotechnology 32:896–902.

    https://doi.org/10.1038/nbt.2931

    • PubMed
    • Google Scholar
88. 1. Robinson MD
    2. McCarthy DJ
    3. Smyth GK

    (2010) edgeR: a bioconductor package for differential expression analysis of digital gene expression data

    Bioinformatics 26:139–140.

    https://doi.org/10.1093/bioinformatics/btp616

    • PubMed
    • Google Scholar
89. 1. Robinson CM
    2. Neary R
    3. Levendale A
    4. Watson CJ
    5. Baugh JA

    (2012) Hypoxia-induced DNA hypermethylation in human pulmonary fibroblasts is associated with Thy-1 promoter methylation and the development of a pro-fibrotic phenotype

    Respiratory Research 13:74.

    https://doi.org/10.1186/1465-9921-13-74

    • PubMed
    • Google Scholar
90. 1. Samanta D
    2. Semenza GL

    (2017) Maintenance of redox homeostasis by hypoxia-inducible factors

    Redox Biology 13:331–335.

    https://doi.org/10.1016/j.redox.2017.05.022

    • PubMed
    • Google Scholar
91. 1. Schödel J
    2. Oikonomopoulos S
    3. Ragoussis J
    4. Pugh CW
    5. Ratcliffe PJ
    6. Mole DR

    (2011) High-resolution genome-wide mapping of HIF-binding sites by ChIP-seq

    Blood 117:e207–e217.

    https://doi.org/10.1182/blood-2010-10-314427

    • PubMed
    • Google Scholar
92. 1. Schunkert H
    2. König IR
    3. Kathiresan S
    4. Reilly MP
    5. Assimes TL
    6. Holm H
    7. Preuss M
    8. Stewart AFR
    9. Barbalic M
    10. Gieger C
    11. Absher D
    12. Aherrahrou Z
    13. Allayee H
    14. Altshuler D
    15. Anand SS
    16. Andersen K
    17. Anderson JL
    18. Ardissino D
    19. Ball SG
    20. Balmforth AJ
    21. Barnes TA
    22. Becker DM
    23. Becker LC
    24. Berger K
    25. Bis JC
    26. Boekholdt SM
    27. Boerwinkle E
    28. Braund PS
    29. Brown MJ
    30. Burnett MS
    31. Buysschaert I
    32. Carlquist JF
    33. Chen L
    34. Cichon S
    35. Codd V
    36. Davies RW
    37. Dedoussis G
    38. Dehghan A
    39. Demissie S
    40. Devaney JM
    41. Diemert P
    42. Do R
    43. Doering A
    44. Eifert S
    45. Mokhtari NEE
    46. Ellis SG
    47. Elosua R
    48. Engert JC
    49. Epstein SE
    50. de Faire U
    51. Fischer M
    52. Folsom AR
    53. Freyer J
    54. Gigante B
    55. Girelli D
    56. Gretarsdottir S
    57. Gudnason V
    58. Gulcher JR
    59. Halperin E
    60. Hammond N
    61. Hazen SL
    62. Hofman A
    63. Horne BD
    64. Illig T
    65. Iribarren C
    66. Jones GT
    67. Jukema JW
    68. Kaiser MA
    69. Kaplan LM
    70. Kastelein JJP
    71. Khaw K-T
    72. Knowles JW
    73. Kolovou G
    74. Kong A
    75. Laaksonen R
    76. Lambrechts D
    77. Leander K
    78. Lettre G
    79. Li M
    80. Lieb W
    81. Loley C
    82. Lotery AJ
    83. Mannucci PM
    84. Maouche S
    85. Martinelli N
    86. McKeown PP
    87. Meisinger C
    88. Meitinger T
    89. Melander O
    90. Merlini PA
    91. Mooser V
    92. Morgan T
    93. Mühleisen TW
    94. Muhlestein JB
    95. Münzel T
    96. Musunuru K
    97. Nahrstaedt J
    98. Nelson CP
    99. Nöthen MM
    100. Olivieri O
    101. Patel RS
    102. Patterson CC
    103. Peters A
    104. Peyvandi F
    105. Qu L
    106. Quyyumi AA
    107. Rader DJ
    108. Rallidis LS
    109. Rice C
    110. Rosendaal FR
    111. Rubin D
    112. Salomaa V
    113. Sampietro ML
    114. Sandhu MS
    115. Schadt E
    116. Schäfer A
    117. Schillert A
    118. Schreiber S
    119. Schrezenmeir J
    120. Schwartz SM
    121. Siscovick DS
    122. Sivananthan M
    123. Sivapalaratnam S
    124. Smith A
    125. Smith TB
    126. Snoep JD
    127. Soranzo N
    128. Spertus JA
    129. Stark K
    130. Stirrups K
    131. Stoll M
    132. Tang WHW
    133. Tennstedt S
    134. Thorgeirsson G
    135. Thorleifsson G
    136. Tomaszewski M
    137. Uitterlinden AG
    138. van Rij AM
    139. Voight BF
    140. Wareham NJ
    141. Wells GA
    142. Wichmann H-E
    143. Wild PS
    144. Willenborg C
    145. Witteman JCM
    146. Wright BJ
    147. Ye S
    148. Zeller T
    149. Ziegler A
    150. Cambien F
    151. Goodall AH
    152. Cupples LA
    153. Quertermous T
    154. März W
    155. Hengstenberg C
    156. Blankenberg S
    157. Ouwehand WH
    158. Hall AS
    159. Deloukas P
    160. Thompson JR
    161. Stefansson K
    162. Roberts R
    163. Thorsteinsdottir U
    164. O'Donnell CJ
    165. McPherson R
    166. Erdmann J

    (2011) Large-scale association analysis identifies 13 new susceptibility loci for coronary artery disease

    Nature Genetics 43:333–338.

    https://doi.org/10.1038/ng.784

    • Google Scholar
93. 1. Shah SH
    2. Roselli A
    3. Honghuang C
    4. Sveinbjornsson L
    5. Fatemifar G

    (2020) Genome-wide association study provides new insights into the genetic architecture and pathogenesis of heart failure

    Nature Communications 11:163.

    https://doi.org/10.1038/s41467-019-13690-5

    • Google Scholar
94. 1. Sigurdsson MI
    2. Saddic L
    3. Heydarpour M
    4. Chang TW
    5. Shekar P
    6. Aranki S
    7. Couper GS
    8. Shernan SK
    9. Muehlschlegel JD
    10. Body SC

    (2017) Post-operative atrial fibrillation examined using whole-genome RNA sequencing in human left atrial tissue

    BMC Medical Genomics 10:25.

    https://doi.org/10.1186/s12920-017-0270-5

    • PubMed
    • Google Scholar
95. 1. Smirnov DA
    2. Brady L
    3. Halasa K
    4. Morley M
    5. Solomon S
    6. Cheung VG

    (2012) Genetic variation in radiation-induced cell death

    Genome Research 22:332–339.

    https://doi.org/10.1101/gr.122044.111

    • PubMed
    • Google Scholar
96. Software

    1. Smith A
    2. Hubley R
    3. Green P

    (2010) RepeatMasker, version 3.0

    Genome Bioinformatics.

    http://www.repeatmasker.org/
97. 1. Smyth GK

    (2004) Linear models and empirical bayes methods for assessing differential expression in microarray experiments

    Statistical Applications in Genetics and Molecular Biology 3:1–25.

    https://doi.org/10.2202/1544-6115.1027

    • PubMed
    • Google Scholar
98. 1. Stephens M

    (2017) False discovery rates: a new deal

    Biostatistics 18:275–294.

    https://doi.org/10.1093/biostatistics/kxw041

    • PubMed
    • Google Scholar
99. 1. Stitziel NO
    2. Stirrups KE
    3. Masca NG
    4. Erdmann J
    5. Ferrario PG
    6. König IR
    7. Weeke PE
    8. Webb TR
    9. Auer PL
    10. Schick UM
    11. Lu Y
    12. Zhang H
    13. Dube MP
    14. Goel A
    15. Farrall M
    16. Peloso GM
    17. Won HH
    18. Do R
    19. van Iperen E
    20. Kanoni S
    21. Kruppa J
    22. Mahajan A
    23. Scott RA
    24. Willenberg C
    25. Braund PS
    26. van Capelleveen JC
    27. Doney AS
    28. Donnelly LA
    29. Asselta R
    30. Merlini PA
    31. Duga S
    32. Marziliano N
    33. Denny JC
    34. Shaffer CM
    35. El-Mokhtari NE
    36. Franke A
    37. Gottesman O
    38. Heilmann S
    39. Hengstenberg C
    40. Hoffman P
    41. Holmen OL
    42. Hveem K
    43. Jansson JH
    44. Jöckel KH
    45. Kessler T
    46. Kriebel J
    47. Laugwitz KL
    48. Marouli E
    49. Martinelli N
    50. McCarthy MI
    51. Van Zuydam NR
    52. Meisinger C
    53. Esko T
    54. Mihailov E
    55. Escher SA
    56. Alver M
    57. Moebus S
    58. Morris AD
    59. Müller-Nurasyid M
    60. Nikpay M
    61. Olivieri O
    62. Lemieux Perreault LP
    63. AlQarawi A
    64. Robertson NR
    65. Akinsanya KO
    66. Reilly DF
    67. Vogt TF
    68. Yin W
    69. Asselbergs FW
    70. Kooperberg C
    71. Jackson RD
    72. Stahl E
    73. Strauch K
    74. Varga TV
    75. Waldenberger M
    76. Zeng L
    77. Kraja AT
    78. Liu C
    79. Ehret GB
    80. Newton-Cheh C
    81. Chasman DI
    82. Chowdhury R
    83. Ferrario M
    84. Ford I
    85. Jukema JW
    86. Kee F
    87. Kuulasmaa K
    88. Nordestgaard BG
    89. Perola M
    90. Saleheen D
    91. Sattar N
    92. Surendran P
    93. Tregouet D
    94. Young R
    95. Howson JM
    96. Butterworth AS
    97. Danesh J
    98. Ardissino D
    99. Bottinger EP
    100. Erbel R
    101. Franks PW
    102. Girelli D
    103. Hall AS
    104. Hovingh GK
    105. Kastrati A
    106. Lieb W
    107. Meitinger T
    108. Kraus WE
    109. Shah SH
    110. McPherson R
    111. Orho-Melander M
    112. Melander O
    113. Metspalu A
    114. Palmer CN
    115. Peters A
    116. Rader D
    117. Reilly MP
    118. Loos RJ
    119. Reiner AP
    120. Roden DM
    121. Tardif JC
    122. Thompson JR
    123. Wareham NJ
    124. Watkins H
    125. Willer CJ
    126. Kathiresan S
    127. Deloukas P
    128. Samani NJ
    129. Schunkert H
    130. Myocardial Infarction Genetics and CARDIoGRAM Exome Consortia Investigators

    (2016) Coding variation in ANGPTL4, LPL, and SVEP1 and the risk of coronary disease

    The New England Journal of Medicine 374:1134–1144.

    https://doi.org/10.1056/NEJMoa1507652

    • PubMed
    • Google Scholar
100. 1. Stone G
     2. Choi A
     3. Meritxell O
     4. Gorham J
     5. Heydarpour M
     6. Seidman CE
     7. Seidman JG
     8. Aranki SF
     9. Body SC
     10. Carey VJ
     11. Raby BA
     12. Stranger BE
     13. Muehlschlegel JD

     (2019) Sex differences in gene expression in response to ischemia in the human left ventricular myocardium

     Human Molecular Genetics 28:1682–1693.

     https://doi.org/10.1093/hmg/ddz014

     • PubMed
     • Google Scholar
101. 1. Strober BJ
     2. Elorbany R
     3. Rhodes K
     4. Krishnan N
     5. Tayeb K
     6. Battle A
     7. Gilad Y

     (2019) Dynamic genetic regulation of gene expression during cellular differentiation

     Science 364:1287–1290.

     https://doi.org/10.1126/science.aaw0040

     • PubMed
     • Google Scholar
102. 1. Unoki M
     2. Nakamura Y

     (2003) Methylation at CpG islands in intron 1 of EGR2 confers enhancer-like activity

     FEBS Letters 554:67–72.

     https://doi.org/10.1016/s0014-5793(03)01092-5

     • PubMed
     • Google Scholar
103. 1. van de Geijn B
     2. McVicker G
     3. Gilad Y
     4. Pritchard JK

     (2015) WASP: allele-specific software for robust molecular quantitative trait locus discovery

     Nature Methods 12:1061–1063.

     https://doi.org/10.1038/nmeth.3582

     • PubMed
     • Google Scholar
104. 1. Wang Y
     2. Ju C
     3. Hu J
     4. Huang K
     5. Yang L

     (2019) PRMT4 overexpression aggravates cardiac remodeling following myocardial infarction by promoting cardiomyocyte apoptosis

     Biochemical and Biophysical Research Communications 520:645–650.

     https://doi.org/10.1016/j.bbrc.2019.10.085

     • PubMed
     • Google Scholar
105. 1. Ward MC
     2. Gilad Y

     (2019) A generally conserved response to hypoxia in iPSC-derived cardiomyocytes from humans and chimpanzees

     eLife 8:e42374.

     https://doi.org/10.7554/eLife.42374

     • PubMed
     • Google Scholar
106. 1. Watson CJ
     2. Collier P
     3. Tea I
     4. Neary R
     5. Watson JA
     6. Robinson C
     7. Phelan D
     8. Ledwidge MT
     9. McDonald KM
     10. McCann A
     11. Sharaf O
     12. Baugh JA

     (2014) Hypoxia-induced epigenetic modifications are associated with cardiac tissue fibrosis and the development of a myofibroblast-like phenotype

     Human Molecular Genetics 23:2176–2188.

     https://doi.org/10.1093/hmg/ddt614

     • PubMed
     • Google Scholar
107. 1. Webb TR
     2. Erdmann J
     3. Stirrups KE
     4. Stitziel NO
     5. Masca NG
     6. Jansen H
     7. Kanoni S
     8. Nelson CP
     9. Ferrario PG
     10. König IR
     11. Eicher JD
     12. Johnson AD
     13. Hamby SE
     14. Betsholtz C
     15. Ruusalepp A
     16. Franzén O
     17. Schadt EE
     18. Björkegren JL
     19. Weeke PE
     20. Auer PL
     21. Schick UM
     22. Lu Y
     23. Zhang H
     24. Dube MP
     25. Goel A
     26. Farrall M
     27. Peloso GM
     28. Won HH
     29. Do R
     30. van Iperen E
     31. Kruppa J
     32. Mahajan A
     33. Scott RA
     34. Willenborg C
     35. Braund PS
     36. van Capelleveen JC
     37. Doney AS
     38. Donnelly LA
     39. Asselta R
     40. Merlini PA
     41. Duga S
     42. Marziliano N
     43. Denny JC
     44. Shaffer C
     45. El-Mokhtari NE
     46. Franke A
     47. Heilmann S
     48. Hengstenberg C
     49. Hoffmann P
     50. Holmen OL
     51. Hveem K
     52. Jansson JH
     53. Jöckel KH
     54. Kessler T
     55. Kriebel J
     56. Laugwitz KL
     57. Marouli E
     58. Martinelli N
     59. McCarthy MI
     60. Van Zuydam NR
     61. Meisinger C
     62. Esko T
     63. Mihailov E
     64. Escher SA
     65. Alver M
     66. Moebus S
     67. Morris AD
     68. Virtamo J
     69. Nikpay M
     70. Olivieri O
     71. Provost S
     72. AlQarawi A
     73. Robertson NR
     74. Akinsansya KO
     75. Reilly DF
     76. Vogt TF
     77. Yin W
     78. Asselbergs FW
     79. Kooperberg C
     80. Jackson RD
     81. Stahl E
     82. Müller-Nurasyid M
     83. Strauch K
     84. Varga TV
     85. Waldenberger M
     86. Zeng L
     87. Chowdhury R
     88. Salomaa V
     89. Ford I
     90. Jukema JW
     91. Amouyel P
     92. Kontto J
     93. Nordestgaard BG
     94. Ferrières J
     95. Saleheen D
     96. Sattar N
     97. Surendran P
     98. Wagner A
     99. Young R
     100. Howson JM
     101. Butterworth AS
     102. Danesh J
     103. Ardissino D
     104. Bottinger EP
     105. Erbel R
     106. Franks PW
     107. Girelli D
     108. Hall AS
     109. Hovingh GK
     110. Kastrati A
     111. Lieb W
     112. Meitinger T
     113. Kraus WE
     114. Shah SH
     115. McPherson R
     116. Orho-Melander M
     117. Melander O
     118. Metspalu A
     119. Palmer CN
     120. Peters A
     121. Rader DJ
     122. Reilly MP
     123. Loos RJ
     124. Reiner AP
     125. Roden DM
     126. Tardif JC
     127. Thompson JR
     128. Wareham NJ
     129. Watkins H
     130. Willer CJ
     131. Samani NJ
     132. Schunkert H
     133. Deloukas P
     134. Kathiresan S
     135. Wellcome Trust Case Control Consortium
     136. MORGAM Investigators
     137. Myocardial Infarction Genetics and CARDIoGRAM Exome Consortia Investigators

     (2017) Systematic evaluation of pleiotropy identifies 6 further loci associated with Coronary Artery Disease

     Journal of the American College of Cardiology 69:823–836.

     https://doi.org/10.1016/j.jacc.2016.11.056

     • PubMed
     • Google Scholar
108. 1. Wei Y
     2. Tenzen T
     3. Ji H

     (2015) Joint analysis of differential gene expression in multiple studies using correlation motifs

     Biostatistics 16:31–46.

     https://doi.org/10.1093/biostatistics/kxu038

     • PubMed
     • Google Scholar
109. 1. Wheeler HE
     2. Shah KP
     3. Brenner J
     4. Garcia T
     5. Aquino-Michaels K
     6. Cox NJ
     7. Nicolae DL
     8. Im HK
     9. GTEx Consortium

     (2016) Survey of the heritability and sparse architecture of gene expression traits across human tissues

     PLOS Genetics 12:e1006423.

     https://doi.org/10.1371/journal.pgen.1006423

     • PubMed
     • Google Scholar
110. Report

     1. WHO

     (2018) World Health Statistics 2018:Monitoring Health for the SDGs

     World Health Organization.

     https://apps.who.int/iris/handle/10665/272596

     • Google Scholar
111. 1. Zhang Y
     2. Liu T
     3. Meyer CA
     4. Eeckhoute J
     5. Johnson DS
     6. Bernstein BE
     7. Nusbaum C
     8. Myers RM
     9. Brown M
     10. Li W
     11. Liu XS

     (2008) Model-based analysis of ChIP-Seq (MACS)

     Genome Biology 9:R137.

     https://doi.org/10.1186/gb-2008-9-9-r137

     • PubMed
     • Google Scholar

Acceptance summary:

This manuscript describes an in vitro system to study inter-individual variation of transcriptional and epigenetic responses to hypoxia using iPSC derived cardiomyocytes. Changes in gene expression are found without changes in chromatin accessibility or methylation, providing insight into key questions about the necessity of epigenetic changes for changes in gene expression. It also identifies regions of the genome responsible for context dependent-changes in gene expression.

Decision letter after peer review:

Thank you for submitting your article "Dynamic effects of genetic variation on gene expression revealed following hypoxic stress in cardiomyocytes" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Patricia Wittkopp as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Senior Editor has drafted this decision to help you prepare a revised submission.

Both reviewers agree that this manuscript is in principle appropriate for publication in eLife without any further experimentation; however, they both also raise a number of concerns about the data analysis and presentation that would need to be addressed satisfactorily before we can consider this work further. Although the typical practice for eLife is to provide authors with a single consolidated review, I felt it would be more beneficial in this case to provide the two reviews as is. They are generally consistent and both explain why they recommend the changes they do in detail.

Reviewer #2:

1) The authors present the range of cardiomyocyte purity. However, I think it would be helpful for this data to be presented in a table. Further, how do the purities correlate between the individual at different time points collected? Could some of the eQTLs identified really be a result of the variation in purity of the differentiations? Investigation and results regarding this possibility would be helpful if included in the manuscript.

2) The inferences made in this paper about the identification of dynamic eQTLs (and their absence) for a given loci rests on the assumption that the power to detect an eQTL is the same in every cell-type by condition. This is addressed in subsection “Dynamic eQTLs are revealed following hypoxia”, where it is clearly show (and acknowledged) that this is not the case. To avoid this problem the authors chose a lenient secondary cut-off at the point where p-values deviate from uniformity. Firstly, it would be helpful to formally test for deviation of uniformity here (rather than appearance on a graph). More importantly, utilising the same p-value threshold when the distribution of p-values is different will lead to variation in the false positive / false negative rates between conditions. This would lead to bias in the identification of dynamic eQTLs, as variation in eQTL detection is a function of the FP/FN rates. Can this be addressed, and it clearly demonstrated that the comparison between conditions is based on the same ability to detect a true eQTL?

3) In subsection “Dynamic eQTLs are revealed following hypoxia”, the authors state that the genes that were not eQTLs in any condition were uniform. However, the figures suggest significant skewing to more significant p-values. Please test for the uniform distribution with a chi-squared test. The authors suggest that they would expect the eQTLs identified in one condition to be uniform in another condition. This relies a the huge assumption that the largest contributing factor to significant eQTLs in each condition is the condition itself, not the cell type and that there will not be a significant number of eQTLs that are consistent across all conditions. However, the authors have provided no evidence to support this claim or this assumption. Please provide the requisite evidence to support these claims and assumptions.

4) Related to this, in paragraph three, it is indicated that based on the FDR calculations, half of the dynamic eQTLs that they identified are not truly dynamic. However, given this is the main conclusion of the paper, it really would be beneficial if some additional work could be done to ascertain the true from false dynamic eQTLs. Particularly if this work is going to be used for any translational research. One option would be to create a null distribution of the effect sizes for genes in each condition and test for enrichment of the "dynamic eQTLs" across the different conditions to test demonstrate whether they are actually identifying truly dynamic eQTLs. Further, I think different FDR cutoffs would have been more appropriate for identifying truly dynamic eQTLs since that is the main focus of the manuscript.

5) The analysis of the GTEx and GWAS overlap is unclear and would benefit from some more careful explanations. Reading this section, I had the following concerns:

i) How were the 14 tissue types randomly selected? And why not include all GTEx tissue types?

ii) How was overlap between dynamic eQTL and GTEx eQTL determined? Many eGenes are regulated by different eQTL – was the overlap based ona. Rank test of p-values between datasets for example? That would help provide confidence it is a true shared effect rather than unlinked eQTL.

iii) Are the allelic effects in the same direction?

iv) Why perform two stage comparison between eQTL and GWAS variants from all traits, then cardiac traits? Was multiple testing corrected for?

v) How do you rule out whether the overlap is just by chance, pleiotropy, linkage, or a causal relationship? These different scenarios have significant impact on the conclusions that are drawn. As it currently stands, it is implied that the overlap is a shared causal relationship.

6) One of the premises of this work is the statement made at the end of paragraph one of the Introduction, where the authors conclude that because the majority GTEx eQTLs are shared across tissues, they probably do not contribute to disease in tissue-specific manners. However, this statement makes two significant assumptions that are not addressed. 1) They assume that the downstream effects and signalling are the same regardless of cell type given they have the same eQTL. However, because different genes and pathways are activated in different tissues, this assumption may not hold true across all tissues. 2) They ignore the ascertainment bias associated with large additive cis-eQTL (which are more likely to overlap between tissues). The mean proportion of expression h2 that can be found to overlap between tissues is 0.02 (2%) in GTEX. The assumption that the remaining 98% are shared between tissues isn't supported by evidence from trans-eQTL mapping or single cell eQTL analyses. In my opinion, clarifying this statement in the Introduction as justification for the work would be beneficial.

7) I found that manuscript in general to be difficult to read and easily interpret the analysis conducted, and conclusions based on specific results. I would suggest a careful eye on editing a review would be very helpful. But here are a few points that I noted specifically

i) The notation used for supplementary figures and tables presented as "Figure X – Supplement Y". It would help to refer to the supplementary figures as "Supplementary Figure X" instead.

ii) Subsection “Dynamic eGenes are enriched for response genes and transcription factors” – it is not always entirely clear whether the authors are referring to transcription factor genes or transcription factors regulating the genes. You could include the word "gene" following "transcription factor" throughout the paragraph.

iii) Defining clearly what response genes and why they are important.

iv) Subsection “Potential limitations of our model” states that the power analysis indicates a 6% power to detect 38% heritability. Firstly, I assume, but it is not clear, that the test is against H0: h2=0? I understand the maths behind this, but that seems a really odd way of explaining power. Why not give a power curve?

v) How are the genomic features tested for enrichment of DARs?

vi) Subsection “Genomic features associated with differentially accessible regions (DARs)” paragraph two and Discussion paragraph four are confusing and unclear. It is hard to understand what the has been done.

vii) Throughout the manuscript, it is difficult to identify which timepoint were used for each assay. Please clarify this when introducing a new assay. Altering Figure 1 to make this clear graphically would help as well.

8) Some of the RIN numbers in RIN table are very low (i.e. 0 for 18852 in condition D). It's not clear if these samples were included in the analyses, but if they were some justification is needed. I worry about impact on the analysis, as essentially it would induce bias in those samples.

9) The section "DNA methylation levels are largely stable following hypoxia and re-oxygenation" is confusing and difficult to read. It is first stated that no DMCpGs were identified and then later stated that they were identified.

10) In subsection “Chromatin accessibility changes following hypoxia and re-oxygenation”, it is stated that ATAC-seq reads were selected that did not show allelic mapping biases. However, later testing for differential accessible regions. This filtering step would remove the ability to effectively identify differentially accessible regions with this step. Could you please either explain the reasoning for removing allelic mapping biases at this step, and/or explain why it would not impact differential accessibility analyses downstream.

11) It is concluded that the data indicate that the "underlying DNA sequence features can affect their epigenetic profile". However, I do not think that the results demonstrate this. They only demonstrate differential methylation in differentially accessible regions which are themselves not associated with genetic variants. I would suggest altering this sentence by removing "that the underlying DNA sequence features of these regions can affect their epigenetic profile"

Reviewer #3:

The authors established an in vitro system to study inter-individual variation of transcriptional and epigenetic responses to hypoxia using iPSC derived cardiomyocytes. Despite the modest sample size, they discovered more than a thousand eQTLs in total, of which 367 are dynamic eQTLs modulating the transcriptional response to hypoxia with implications for disease susceptibility.

Although the analyses demonstrated in the manuscript are solid and the results are reliable, authors failed to present their results effectively. I would recommend the authors to rewrite the paper to increase the accuracy and readability in general (as I suggested below).

In addition, I have one particular concern regarding the following important (but vague) claim presented in this paper, that is the transcriptional response to hypoxia is not mediated through epigenetic changes. Although the authors repeated the claim multiple times throughout the paper, it is hugely confounded by the power to detect epigenetic variations between conditions/individuals (i.e., DARs, caQTLs and DMCpGs). If the authors would keep the claim, they need to address the power issue appropriately.

Here are my specific comments:

1) Figure 1A should explain the overview of the experiment in more detail. There are three replicates for the three of 15 donors (mentioned in the main text) which are hidden in the panel. It is also nice to combine panel A and D in one panel to illustrate the whole experiment. According to the main text, one donor was dropped in ATAC-seq and two donors were dropped in methylation data which should also be shown in the panel.

2) It would be nice to summarise the 4 x 21 RNA-seq samples in the main figure. The authors used a linear mixed model for differential expression (DE) analysis. I'm wondering if the model can be used to perform a variance components analysis to decompose the total variance (for each gene) into “condition”, “individual”, “replication”, and residual variances (where the condition and replication are also treated as random effects). The estimated variance parameters across all genes are then summarised by a boxplot or similar where the residual variance is standardised as 1.0. The authors could repeat the analysis for ATAC-seq and methylation data to compare the difference of explained variances between molecular phenotypes. The analysis will provide clear interpretation on why detectable DARs, caQTLs and DMCpGs are so few. The similar analysis was performed e.g., in (Kilpinen et al., 2017).

3) Figure 2H is not so informative. Instead, the authors could extend the heatmap in Figure 2F to show how much of eGenes are classified into baseline eGenes, dynamic eGenes and shared eGenes by showing all genes with eQTLs at least in one condition. The authors could also add additional columns next to the heatmap showing which eGenes are response genes and/or transcription factors to support the enrichment analysis in the section “Dynamic eGenes are enriched for response genes and…”.

4) The DE analysis result is frequently used in the manuscript, but there is no summary data shown in the main text (e.g., volcano plots, DE genes shared/not shared with different conditions, etc.). It would be ideal to have an independent figure for DE analysis. Currently, Figure 2A-C are the result of DE analysis which are inconsistent with the figure title (i.e., “Hundreds of dynamic eQTLs are revealed…”).

5) How many dynamic eGenes are overlapping with the baseline eGenes? The Chi-square test performed in the section “Dynamic eGenes are enriched for response genes and transcription factors” seems inappropriate because the authors double-count the same genes in baseline eGenes as the background distribution. They need to subtract dynamic eGenes from the baseline eGenes to perform a proper Chi-square test.

6) The authors performed the enrichment analysis of transcription factors (TFs) in dynamic eGenes which could be confounded by response genes, because dynamic eGenes are likely to be response genes. In addition, even if dynamic eGenes are enriched with TFs, what is the biological interpretation? Why do regulatory variants affect TF dynamic eGenes more than non-TF dynamic eGenes?

7) It is often the case that the main figures are used to show anecdotal examples of biological statements raised in the main text and important data/results are hidden in the supplementary figures. The main figures must be used to support the author's claims in general (without cherry-picking an example gene). For example, the authors stated “most DARs have returned to baseline levels of accessibility by the first re-oxygenation condition.” with just numbers of DARs and the FOXO1 example which do not strongly support the claim. The authors should perform and highlight global analyses to support each claim from different angles. In fact, the claim was partly supported by the supplementary figure 3 (supplement 4A-C). Therefore, the authors first consider rearranging the figures to better support each claim made in the main text. Then, they can clarify which analysis is needed to make a strong statement.

8) The statement “These results suggests that gene expression changes…” was supported only by a limited number of DARs and not totally addressed yet. Recently, the “gene activity score” was proposed in the single-cell ATAC-seq analysis to overcome the sparseness of chromatin data (Granja et al., 2020). I'm wondering if it improves the power to detect DARs in the paper. This quantification also allows the authors to easily colocalise DARs with eGenes/response genes to prove (or disprove) the statement.

9) It is very surprising that CHT found only tens of caQTLs, while it mapped more than a thousand eQTLs. The authors should comment on that. If the authors use only variants overlapping with the accessibility region, does the power to detect caQTL increase? It seems 25Kb cis-region is too broad.

10) The statement “Changes in DNA methylation are therefore…” is hugely confounded by the power to detect DMCpGs with the higher noise rate in DNA methylation arrays. It is not conclusive with the modest sample size.

11) The authors should perform a colocalisation analysis (Giambartolomei et al., 2014) if GWAS summary statistics are available. If not, the authors are still able to calculate linkage disequilibrium (LD) between the GWAS lead variant and eQTL lead variant from CHT to make sure the GWAS locus and the eQTL are reasonably colocalised. The authors should show locuszoom plots in conjunction with boxplots in Figure 6 to demonstrate both the magnitude of statistical significance of eQTLs and LD (r² index) with GWAS lead variants.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for submitting your article "Dynamic effects of genetic variation on gene expression revealed following hypoxic stress in cardiomyocytes" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Patricia Wittkopp as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.

As stated by the reviewers in response to the first submission, the paper presents a timely, interesting and relevant advance. The initial comments raised have been largely addressed, however both reviewers continue raise concerns regarding the eQTL and GWAS overlap / colocalization analysis. We agree that these concerns are fundamental and should be addressed.

First, concerning the overlap with GTEX (R2), it seems entirely doable to perform this analysis on the level of individual variants, which also permit assessing consistency of effect directions, etc.

Second, both reviewers comment on the GWAS overlap/colocalization analysis. While we acknowledge that the resolution to definitely identify colocalization events will be limited, the authors should extend their analysis and/or town down the conclusions. The approach proposed by reviewer 3 seems sensible and viable, at least for a selected set of loci.

For reference, we have also included the complete reviewer reports.

Reviewer #2:

The authors have addressed most of the comments raised, although there are a couple of points that I feel would be beneficial to investigate further.

1) GTEX overlap – in the original review, a suggestion of testing for concordance of allelic effect direction be shown for the replication of the eQTL identified in this work and GTEX tissues. It is unclear what “analysis at the level of the gene” exactly means, but given that many genes are controlled by independent eQTL (loci) in difference tissues/cell types, overlapping an observation of eGenes doesn't allow conclusions of shared eQTL to be determined. Testing for the concordance of allelic effects for eSNPs identified from your analysis with GTEX eSNP results would provide much stronger evidence to support this conclusion.

2) GWAS overlap – "We next chose to analyze the data using a less stringent gene-based approach using GWAS traits of interest. The second analysis allows for the identification of potentially interesting regions where the causal SNP may be different in the eQTL and GWAS data" I am not sure what this means biologically and what conclusions can be drawn from this. What does it mean if there is no overlap in the SNPs between a eQTL and GWAS loci but that they are both in the same gene region? Is there an enrichment over what would be expected by chance? the statements "We did not mean to imply that the overlap is a shared causal relationship. We have changed the language to reflect this. Our goal is to illustrate that the dynamic eQTLs we identified may play a role in higher-level phenotypes." These sentences seems in opposition to one another.

Reviewer #3:

The authors have now addressed most of my questions and concerns. Yet, the result of GWAS colocalisation is still in question. I have previously suggested to perform a proper colocalisation analysis using COLOC package or to report the linkage disequilibrium index (r² value) between GWAS index variant and eQTL lead variant. Unfortunately the authors just provided only P-values for associations which does not support the GWAS locus and eQTL are colocalised. The authors have to expand the cis-window (up to 1Mb) for a handful number of genes overlapping with GWAS loci to perform a proper colocalisation analsysis. It is crucial to provide the posterior probability of colocalisation for such a small number of GWAS overlapping genes.

Reviewer #2:

1) The authors present the range of cardiomyocyte purity. However, I think it would be helpful for this data to be presented in a table. Further, how do the purities correlate between the individual at different time points collected? Could some of the eQTLs identified really be a result of the variation in purity of the differentiations? Investigation and results regarding this possibility would be helpful if included in the manuscript.

Please see Supplementary file 1 as indicated in the text for the purity of each cell line.

We did not measure purity at different time points, only at day 30 +/- 1 day, before the treatment on day 31-32. Because we split treatment cultures from a single differentiation culture, and all experiments were concluded within 48h after purity was determined by flow cytometry, it is unlikely that there are purity differences between conditions.

We have made this clearer in the Results and we have added a schematic of the experimental workflow as Figure 1—figure supplement 1. We have also amended Figure 2—figure supplement 4 to include cardiomyocyte marker gene expression levels across all conditions instead of just condition A as a representative condition. This data illustrates that there is no systematic bias in cardiomyocyte marker gene expression levels across conditions suggesting that the purity is indeed similar across conditions.

2) The inferences made in this paper about the identification of dynamic eQTLs (and their absence) for a given loci rests on the assumption that the power to detect an eQTL is the same in every cell-type by condition. This is addressed in subsection “Dynamic eQTLs are revealed following hypoxia”, where it is clearly show (and acknowledged) that this is not the case. To avoid this problem the authors chose a lenient secondary cut-off at the point where p-values deviate from uniformity. Firstly, it would be helpful to formally test for deviation of uniformity here (rather than appearance on a graph). More importantly, utilising the same p-value threshold when the distribution of p-values is different will lead to variation in the false positive / false negative rates between conditions. This would lead to bias in the identification of dynamic eQTLs, as variation in eQTL detection is a function of the FP/FN rates. Can this be addressed, and it clearly demonstrated that the comparison between conditions is based on the same ability to detect a true eQTL?

We tested for deviation from uniform distributions and added that data to the figure legend of Figure 3—figure supplement 2A. We believe that the assumption that power to detect eQTLs is similar in all conditions is a reasonable one as we now note because: (i) the cells originated from the same differentiation culture for each individual so the only factor distinguishing these samples should be the condition. (ii) We have the same number of individuals in each condition. (iii) RNA was extracted from all time points of the same individual at the same time (Supplementary file 1). (iv) RIN scores are similar across conditions (Supplementary file 1). (v) Sequencing read counts are similar across conditions (Figure 2—figure supplement 1). Moreover, the distribution of eQTL effect sizes across conditions is similar (new supplementary figure: Figure 3—figure supplement 1), and the values are correlated, further suggesting we have similar power to detect eQTLs across conditions.

The reviewer is correct that reliance on an arbitrary P-value when the distributions are different will lead to a different FDR. For that reason, our first – more stringent – cutoff is based on an FDR value, not a P-value. The secondary cutoff is based on a nominal P-value, which was determined based on visualization of the distributions. Importantly, we provide the FDR associated with different secondary P-value cutoffs in reviewer comment four below, thereby addressing the reviewer’s concern.

3) In subsection “Dynamic eQTLs are revealed following hypoxia”, the authors state that the genes that were not eQTLs in any condition were uniform. However, the figures suggest significant skewing to more significant p-values. Please test for the uniform distribution with a chi-squared test. The authors suggest that they would expect the eQTLs identified in one condition to be uniform in another condition. This relies a the huge assumption that the largest contributing factor to significant eQTLs in each condition is the condition itself, not the cell type and that there will not be a significant number of eQTLs that are consistent across all conditions. However, the authors have provided no evidence to support this claim or this assumption. Please provide the requisite evidence to support these claims and assumptions.

The experiments were concluded within 48 hours of testing a representative sample from each individual for purity. It does not seem reasonable that during that time the composition of the cultures has changed, certainly not in a way that will create a biased cell composition association with a specific treatment condition. We now acknowledge in the text that our conclusion relies on this assumption, and explain why we believe that it is quite a reasonable assumption to make.

We have tested for deviation from a uniform distribution and indicated the statistics in the figure legend of Figure 3—figure supplement 2B, and in the main text. However, given that the outcome of a chi-square test (rejection/non-rejection) could depend on how the data are binned, and it is not obvious how to do so in this case, we used the Kolmogorov–Smirnov test to determine whether our P-value distribution deviates from a uniform distribution as a simpler alternative.

4) Related to this, in paragraph three, it is indicated that based on the FDR calculations, half of the dynamic eQTLs that they identified are not truly dynamic. However, given this is the main conclusion of the paper, it really would be beneficial if some additional work could be done to ascertain the true from false dynamic eQTLs. Particularly if this work is going to be used for any translational research. One option would be to create a null distribution of the effect sizes for genes in each condition and test for enrichment of the "dynamic eQTLs" across the different conditions to test demonstrate whether they are actually identifying truly dynamic eQTLs. Further, I think different FDR cutoffs would have been more appropriate for identifying truly dynamic eQTLs since that is the main focus of the manuscript.

As the reviewer appreciates, it is not obvious how to distinguish between the true and false positives in this case. We supposed that one can use priors based on gene functions, but this approach seems prone to errors. We have now determined the number of dynamic eQTLs at different statistical cutoffs and their associated FDR (P=0.5: 117 dynamic eQTLs with 30% FDR; P=0.75: 53 eQTLs with 13% FDR; P=0.8: 37 eQTLs with 11% FDR;, P=0.9: 19 eQTLs with 7% FDR). We have added the results of one of these more stringent P-value cutoffs, where the FDR is quite low (7%) to the main text. We also now provide a supplementary figure highlighting six additional examples of dynamic eQTLs (Figure 3—figure supplement 3) indicating that our approach identifies dynamic eQTLs as expected.

5) The analysis of the GTEx and GWAS overlap is unclear and would benefit from some more careful explanations. Reading this section, I had the following concerns:

i) How were the 14 tissue types randomly selected? And why not include all GTEx tissue types?

This number and the specific tissues were chosen in an arbitrary way. We were looking for a modest number of diverse tissue types to interrogate. This analysis illustrates that 98% of dynamic eGenes are eGenes in at least one of 14 tissues suggesting we are reaching saturation already. Adding more tissues therefore can’t add much more insight because we are already at 98% overlap. We clarified that the choice of tissues was arbitrary – effectively random in the text.

ii) How was overlap between dynamic eQTL and GTEx eQTL determined? Many eGenes are regulated by different eQTL – was the overlap based ona. Rank test of p-values between datasets for example? That would help provide confidence it is a true shared effect rather than unlinked eQTL.

As indicated in the Materials and methods this analysis was performed at a gene level; i.e. a gene’s expression level is known to be affected by genetic variation (eGene) in GTEx data. We asked how often do these GTEx eGenes overlap with our eGenes. We have made the approach we chose more explicit in the Results.

A more nuanced SNP-based analysis is actually quite challenging to perform and will require a different approach. As we are uncertain what the potential insight would be from such an analysis, on balance, we decided against it.

iii) Are the allelic effects in the same direction?

As the analysis we performed was based on a gene level we are unable to compare the genotype/allelic effects. There is overwhelming overlap between our findings and GTEx. While we agree that a more nuanced approach will provide more details, we argue that this is a minor point in the paper, and further analysis in this section will not substantially add insight.

iv) Why perform two stage comparison between eQTL and GWAS variants from all traits, then cardiac traits? Was multiple testing corrected for?

To identify genomic loci that may be relevant for complex phenotypes, we performed two independent analyses to intersect our eQTL data with data from GWAS studies. We chose two different analytical approaches to identify as many potentially interesting loci as possible. We first chose to analyze the data using an unbiased, stringent SNP-based approach using all available GWAS traits. We next chose to analyze the data using a less stringent gene-based approach using GWAS traits of interest. The second analysis allows for the identification of potentially interesting regions where the causal SNP may be different in the eQTL and GWAS data.

In both of these analyses we searched for direct overlaps of either SNPs or genes. We did not correct for multiple testing. In line with the comment below we have changed the language to more accurately reflect the analysis we performed.

We have also now included a locus-specific co-localization-based analysis as requested by reviewer three. We have clarified the choice of these three different approaches in the text.

v) How do you rule out whether the overlap is just by chance, pleiotropy, linkage, or a causal relationship? These different scenarios have significant impact on the conclusions that are drawn. As it currently stands, it is implied that the overlap is a shared causal relationship.

We did not mean to imply that the overlap is a shared causal relationship. We have changed the language to reflect this. Our goal is to illustrate that the dynamic eQTLs we identified may play a role in higher-level phenotypes.

6) One of the premises of this work is the statement made at the end of paragraph one of the Introduction, where the authors conclude that because the majority GTEx eQTLs are shared across tissues, they probably do not contribute to disease in tissue-specific manners. However, this statement makes two significant assumptions that are not addressed. 1) They assume that the downstream effects and signalling are the same regardless of cell type given they have the same eQTL. However, because different genes and pathways are activated in different tissues, this assumption may not hold true across all tissues. 2) They ignore the ascertainment bias associated with large additive cis-eQTL (which are more likely to overlap between tissues). The mean proportion of expression h2 that can be found to overlap between tissues is 0.02 (2%) in GTEX. The assumption that the remaining 98% are shared between tissues isn't supported by evidence from trans-eQTL mapping or single cell eQTL analyses. In my opinion, clarifying this statement in the Introduction as justification for the work would be beneficial.

As the reviewer notes this is a complex issue, which was overly simplified in our statement about the consequences of sharing of eQTLs across tissues. We have removed this sentence and instead include a statement about how a variety of approaches are necessary to identify the molecular basis of disease-associated loci.

7) I found that manuscript in general to be difficult to read and easily interpret the analysis conducted, and conclusions based on specific results. I would suggest a careful eye on editing a review would be very helpful. But here are a few points that I noted specifically

We have addressed each point below.

i) The notation used for supplementary figures and tables presented as "Figure X – Supplement Y". It would help to refer to the supplementary figures as "Supplementary Figure X" instead.

The figures have been referenced in the format specified by eLife.

ii) Subsection “Dynamic eGenes are enriched for response genes and transcription factors” – it is not always entirely clear whether the authors are referring to transcription factor genes or transcription factors regulating the genes. You could include the word "gene" following "transcription factor" throughout the paragraph.

By “transcription factor” we mean any gene that is annotated as a transcription factor by Lambert et al., not transcription factors that regulate a particular set of genes. This is indicated in the Materials and methods but we have made this more explicit in the Results.

iii) Defining clearly what response genes and why they are important.

We should have explicitly defined response genes. We have now done so in the text.

iv) Subsection “Potential limitations of our model” states that the power analysis indicates a 6% power to detect 38% heritability. Firstly, I assume, but it is not clear, that the test is against H0: h2=0? I understand the maths behind this, but that seems a really odd way of explaining power. Why not give a power curve?

We clarify that the test is against H₀: λ = 0, where λ is the standardized effect size.

This corresponds to the typical regression-based approach to map QTLs, assuming that genotypes and phenotypes have been standardized.

We plotted the power curve as a function of standardized effect size, fixing the sample size to the size of the current study, in Figure 7—figure supplement 3 (green curve).

Assuming a single causal variant, the standardized effect size be used to derive the heritability explained by that variant.

We reported power to detect an effect of a given true effect size, holding the sample size fixed to the size of the current study. We gave the true effect size in terms of heritability explained, rather than giving the standardized effect size, in order to make the results of the power analysis easier to interpret. The standardized effect size does not correspond to log fold change, which would be an alternative, interpretable choice of units.

In particular, the median eQTL heritability over all genes is h² = 0.16, as previously reported by Gusev et al., 2016 and Wheeler et al., 2016. This reference point makes clear the relative power of this study compared to the GTEx and DGN studies for eQTL mapping.

We have added the Gusev and Wheeler studies for reference to the Discussion.

v) How are the genomic features tested for enrichment of DARs?

We have added the statistical test used to identify promoter- and enhancer-associated marks (chi-squared) to the Results. We have also added more details in the Materials and methods to describe how MEME identifies enriched motifs.

vi) Subsection “Genomic features associated with differentially accessible regions (DARs)” paragraph two and Discussion paragraph four are confusing and unclear. It is hard to understand what the has been done.

We have re-written these sections to make them clearer.

vii) Throughout the manuscript, it is difficult to identify which timepoint were used for each assay. Please clarify this when introducing a new assay. Altering Figure 1 to make this clear graphically would help as well.

We used all time points in all assays and differences between all timepoints were tested in all assays.

We have clarified this in the text, re-designed Figures 1 and 2 and added a supplementary figure with the experimental setup (Figure 1—figure supplement 1) to illustrate this point.

8) Some of the RIN numbers in RIN table are very low (i.e. 0 for 18852 in condition D). It's not clear if these samples were included in the analyses, but if they were some justification is needed. I worry about impact on the analysis, as essentially it would induce bias in those samples.

The median RIN score across 84 samples is 9.1 (9.4 for A, 8.7 for B, 9.1 for C and 9.2 for D). All RIN scores are >8, except 18852A which is 6 and 18852D which is listed as 0. The RIN score for 18852D could not be calculated by the bioanalyzer software, despite the two ribosomal subunits evident in the sample. We should have listed this score as “not detected” not “0”. Given that we have multiple replicates for 18852, and that only one replicate can be used in the eQTL and differential expression analysis, we chose not to include this replicate in downstream analysis.

We have added clarification text to where we mention RIN scores in the Materials and methods, changed the RIN score of 18852D to “not detected” in the table to more accurately reflect this score, and plotted the RIN values in Figure 2—figure supplement 1A to highlight this data.

9) The section "DNA methylation levels are largely stable following hypoxia and re-oxygenation" is confusing and difficult to read. It is first stated that no DMCpGs were identified and then later stated that they were identified.

We did not find any differentially methylated CpGs when considering all CpGs. We did find a handful of differentially methylated CpGs when restricting our analysis to a smaller set of CpGs (those in CpG islands etc).

We have made this distinction, and the rationale behind it, clearer in the text. We also now do not use the term DMCpGs when considering all CpGs. We only introduce the term when we subset CpGs and actually identify differences between conditions. This distinction is also now illustrated in re-designed Figure 6.

10) In subsection “Chromatin accessibility changes following hypoxia and re-oxygenation”, it is stated that ATAC-seq reads were selected that did not show allelic mapping biases. However, later testing for differential accessible regions. This filtering step would remove the ability to effectively identify differentially accessible regions with this step. Could you please either explain the reasoning for removing allelic mapping biases at this step, and/or explain why it would not impact differential accessibility analyses downstream.

We treated the ATAC-seq data as we did the RNA-seq data i.e. used the WASP algorithm to reprocess mapped reads to reduce reference bias as described in the Materials and methods. This step does not impact the ability to identify differentially accessible regions across conditions, it merely removes ambiguous reads. We have clarified in the Results that we only maintain ATAC-seq reads that map unambiguously to the genome.

11) It is concluded that the data indicate that the "underlying DNA sequence features can affect their epigenetic profile". However, I do not think that the results demonstrate this. They only demonstrate differential methylation in differentially accessible regions which are themselves not associated with genetic variants. I would suggest altering this sentence by removing "that the underlying DNA sequence features of these regions can affect their epigenetic profile"

We have altered this sentence as suggested.

Reviewer #3:

The authors established an in vitro system to study inter-individual variation of transcriptional and epigenetic responses to hypoxia using iPSC derived cardiomyocytes. Despite the modest sample size, they discovered more than a thousand eQTLs in total, of which 367 are dynamic eQTLs modulating the transcriptional response to hypoxia with implications for disease susceptibility.

Although the analyses demonstrated in the manuscript are solid and the results are reliable, authors failed to present their results effectively. I would recommend the authors to rewrite the paper to increase the accuracy and readability in general (as I suggested below).

In addition, I have one particular concern regarding the following important (but vague) claim presented in this paper, that is the transcriptional response to hypoxia is not mediated through epigenetic changes. Although the authors repeated the claim multiple times throughout the paper, it is hugely confounded by the power to detect epigenetic variations between conditions/individuals (i.e., DARs, caQTLs and DMCpGs). If the authors would keep the claim, they need to address the power issue appropriately.

We used 13-14 biological replicates (different individuals) in our study. Previous studies have generally used just 1-6 biological replicates to identify differences in chromatin profiles or DNA methylation levels under conditions of stress (Aref-Eshghi, Am J Physio, Cell Physiol, 2020; Marr, Plos Path., 2014; Pacis, GR, 2015). It is also our experience that when we measure both gene expression and DNA methylation changes between sample groups, the number of DNA methylation changes often far outnumber the gene expression changes (Burrows/Banovich, Plos Genetics, 2016; Natri/Bobowik, Plos Genetics, 2020). Lack of power is always a reasonable explanation, but compared to previous studies, our observations are quite surprising. Nevertheless, we have now softened our claims about the lack of differences we observed in the Results section, and removed some of the conclusions we drew.

Here are my specific comments:

1) Figure 1A should explain the overview of the experiment in more detail. There are three replicates for the three of 15 donors (mentioned in the main text) which are hidden in the panel. It is also nice to combine panel A and D in one panel to illustrate the whole experiment. According to the main text, one donor was dropped in ATAC-seq and two donors were dropped in methylation data which should also be shown in the panel.

We have re-designed Figure 1 to combine A and D and now indicate exactly which samples were collected.

2) It would be nice to summarise the 4 x 21 RNA-seq samples in the main figure. The authors used a linear mixed model for differential expression (DE) analysis. I'm wondering if the model can be used to perform a variance components analysis to decompose the total variance (for each gene) into “condition”, “individual”, “replication”, and residual variances (where the condition and replication are also treated as random effects). The estimated variance parameters across all genes are then summarised by a boxplot or similar where the residual variance is standardised as 1.0. The authors could repeat the analysis for ATAC-seq and methylation data to compare the difference of explained variances between molecular phenotypes. The analysis will provide clear interpretation on why detectable DARs, caQTLs and DMCpGs are so few. The similar analysis was performed e.g., in (Kilpinen et al., 2017).

We have now created an extra main text figure where we show the differential expression analysis (Figure 2). These results are largely in line with the findings from our previous work (Ward and Gilad, 2019), which is why we had not initially highlighted them in this paper.

Using a linear mixed model, we calculated the proportion of variance explained by individual, condition, replicate (we only have replicate data for the RNA-seq and DNA methylation data), and the residual variance. A boxplot of variance components for all three phenotypes is now included as Figure 6—figure supplement 4. We now note in the main text that we observe a lower proportion of variance explained by “condition” in the DNA methylation data (0.8%) and chromatin accessibility data (3%) compared to the gene expression data (6%, t-test, P < 2x10^-16), and a corresponding increase in the residual variance. Comparing the gene expression data with the DNA methylation data in particular, reveals that while there is approximately an order of magnitude lower variance in the “condition” component in the DNA methylation data, there is a similar proportion of variance explained by “individual” (44% for gene expression, 40% for DNA methylation, and 28% for chromatin accessibility). These results suggest that noise alone cannot explain the relatively small number of chromatin accessibility and DNA methylation changes.

In terms of the low number of caQTLs, we specified in the Results section that we restricted tested SNPs to those within the 128,672 open chromatin regions we identified. We chose to do this so that any effects would be interpretable i.e. a SNP affects a transcription factor binding site, altering accessibility of the region. This means we were testing many fewer variants than the eQTL analyses. We may well be underpowered to detect these effects in our modest sample size. We have added text to the manuscript to make this point more explicit.

3) Figure 2H is not so informative. Instead, the authors could extend the heatmap in Figure 2F to show how much of eGenes are classified into baseline eGenes, dynamic eGenes and shared eGenes by showing all genes with eQTLs at least in one condition. The authors could also add additional columns next to the heatmap showing which eGenes are response genes and/or transcription factors to support the enrichment analysis in the section “Dynamic eGenes are enriched for response genes and…”.

We appreciate the suggestion on how to make this figure more informative. Unfortunately this didn’t work in practice partly because of the distinction between eSNP-based and eGene-based analyses and partly because of the available software. Instead, to address the comment we have now made a clear distinction between SNP-based results and eGene-based results, added an extra panel with a heatmap which includes all eQTLs identified in all conditions to provide a more global view of the data and to show how much “sharing” there is between conditions (Figure 3B), and re-designed a panel to include the response gene and TF data (Figure 3F).

4) The DE analysis result is frequently used in the manuscript, but there is no summary data shown in the main text (e.g., volcano plots, DE genes shared/not shared with different conditions, etc.). It would be ideal to have an independent figure for DE analysis. Currently, Figure 2A-C are the result of DE analysis which are inconsistent with the figure title (i.e., “Hundreds of dynamic eQTLs are revealed…”).

We have made a new Figure 2 which displays pairwise and joint tests in the RNA-seq data using panels from the original Figure 2 and supplement.

5) How many dynamic eGenes are overlapping with the baseline eGenes? The Chi-square test performed in the section “Dynamic eGenes are enriched for response genes and transcription factors” seems inappropriate because the authors double-count the same genes in baseline eGenes as the background distribution. They need to subtract dynamic eGenes from the baseline eGenes to perform a proper Chi-square test.

This was an error on our part. As the reviewer noted a subset of dynamic eGenes are present at baseline – these are the 37 suppressed eQTLs. We have now excluded these genes to compare only baseline eGenes and induced eGenes so that this is now a mutually exclusive set. This in fact increases the enrichment of both response genes and transcription factors, which we now report in the text.

6) The authors performed the enrichment analysis of transcription factors (TFs) in dynamic eGenes which could be confounded by response genes, because dynamic eGenes are likely to be response genes. In addition, even if dynamic eGenes are enriched with TFs, what is the biological interpretation? Why do regulatory variants affect TF dynamic eGenes more than non-TF dynamic eGenes?

We have now investigated the enrichment of TFs in induced eGenes that are response genes compared to induced eGenes that are non-response genes. TFs amongst induced eGenes are more likely to be response genes than non-response genes (P = 0.02); however there is no difference in the proportion of response genes between TFs amongst baseline eGenes and TFs amongst induced eGenes (38% vs. 39%). We have added these results to the main text.

The enrichment we observe in TFs amongst our induced eGenes suggests that latent genetic variation can have multiple downstream effects on gene expression including gene targets of TF genes. This could provide a mechanism for the appearance of induced eGenes that are neither response genes nor TF genes. We have added the implication of our findings to the main text.

7) It is often the case that the main figures are used to show anecdotal examples of biological statements raised in the main text and important data/results are hidden in the supplementary figures. The main figures must be used to support the author's claims in general (without cherry-picking an example gene). For example, the authors stated “most DARs have returned to baseline levels of accessibility by the first re-oxygenation condition.” with just numbers of DARs and the FOXO1 example which do not strongly support the claim. The authors should perform and highlight global analyses to support each claim from different angles. In fact, the claim was partly supported by the supplementary figure 3 (supplement 4A-C). Therefore, the authors first consider rearranging the figures to better support each claim made in the main text. Then, they can clarify which analysis is needed to make a strong statement.

We have re-designed Figure 4 and 6 to include the genome-wide analyses performed for the ATAC-seq and DNA methylation data that were previously shown in the supplementary figures, and we have edited the text accordingly.

8) The statement “These results suggests that gene expression changes…” was supported only by a limited number of DARs and not totally addressed yet. Recently, the “gene activity score” was proposed in the single-cell ATAC-seq analysis to overcome the sparseness of chromatin data (Granja et al., 2020). I'm wondering if it improves the power to detect DARs in the paper. This quantification also allows the authors to easily colocalise DARs with eGenes/response genes to prove (or disprove) the statement.

The referenced statement refers to the identification of caQTLs not DARs. However, given the other reviewer comments we have removed this particular statement and softened the language in this section. The referenced gene activity score is for aggregating regulatory regions across single cells. We aggregate open chromatin regions across individuals within each condition (by merging peaks of read enrichment) and then look for differences in read count between conditions and individuals. It is not clear to us how this tool would work with our bulk data, or how it would improve power to detect differences between conditions.

9) It is very surprising that CHT found only tens of caQTLs, while it mapped more than a thousand eQTLs. The authors should comment on that. If the authors use only variants overlapping with the accessibility region, does the power to detect caQTL increase? It seems 25Kb cis-region is too broad.

As we stated in the Results and in the Materials and methods, we tested for caQTLs as we did for eQTLs except we did not use a 25 kb window, we restricted SNPs to those within the 128,672 open chromatin regions.

10) The statement “Changes in DNA methylation are therefore…” is hugely confounded by the power to detect DMCpGs with the higher noise rate in DNA methylation arrays. It is not conclusive with the modest sample size.

We would argue that most studies aimed at identifying differences between two conditions do not have 13 biological replicates. Nevertheless, we have removed this statement.

11) The authors should perform a colocalisation analysis (Giambartolomei et al., 2014) if GWAS summary statistics are available. If not, the authors are still able to calculate linkage disequilibrium (LD) between the GWAS lead variant and eQTL lead variant from CHT to make sure the GWAS locus and the eQTL are reasonably colocalised. The authors should show locuszoom plots in conjunction with boxplots in Figure 6 to demonstrate both the magnitude of statistical significance of eQTLs and LD (r² index) with GWAS lead variants.

We only have a small number of overlapping results here, so we chose an alternative approach: We were able to obtain full summary statistics from the CARDIoGRAMplusC4D Consortium for coronary artery disease (CAD) and myocardial infarction (MI). From both significant and dynamic eQTLs we tested the lead eSNP against the full GWAS results and identified two new associations. The first is an association between a dynamic eQTL (regulating EIF2B2) and CAD. The p-value of the dynamic eSNP is 4.7x10^-5 in the GWAS and the most significant SNP in the GWAS peak has a p-value of 9.9x10^-8. As we tested a relatively narrow window around genes in our eQTL analysis, our summary statistics do not overlap the full peak. We also identified an association between a significant eQTL in condition B that is significantly associated with MI (CARM1). We have added both of these loci to the manuscript, and now highlight EIF2B2 with locus zoom plots in Figure 7 and CARM1 in Figure 7—figure supplement 2.

The two original loci that we highlighted were RNF166, associated with varicose veins, and ZC3HC1, associated with MI. Varicose vein GWAS summary statistics were only available for the significant associations. We include a locus zoom plot of the RNF166 region with only the significant GWAS hits (Author response image 1A). The GWAS region around ZC3HC1 does not show a clear peak despite being a known MI locus (Author response image 1B). We felt that neither of these locus zoom plots are striking as examples to illustrate in the main paper. However, we believe that these are nevertheless potentially two interesting loci worth mentioning. We have moved the original RNF166 and ZC3HC1 loci eQTL boxplots to Figure 7—figure supplement 1.

Author response image 1

Download asset Open asset

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Reviewer #2:

The authors have addressed most of the comments raised, although there are a couple of points that I feel would be beneficial to investigate further.

1) GTEX overlap – in the original review, a suggestion of testing for concordance of allelic effect direction be shown for the replication of the eQTL identified in this work and GTEX tissues. It is unclear what “analysis at the level of the gene” exactly means, but given that many genes are controlled by independent eQTL (loci) in difference tissues/cell types, overlapping an observation of eGenes doesn't allow conclusions of shared eQTL to be determined. Testing for the concordance of allelic effects for eSNPs identified from your analysis with GTEX eSNP results would provide much stronger evidence to support this conclusion.

We have now performed the GTEx overlap analysis of dynamic and shared eQTLs at the level of eSNP-eGene pairs rather than eGenes. 326 of 367 dynamic eSNPs are tested in GTEx. 32 dynamic eSNP-eGene pairs are eQTLs in left ventricle heart tissue (9.8%), and 140 are eQTLs in at least one of the other 49 tissues assayed by GTEx (42.9%). For the shared eQTL analysis we selected the 61 eGenes present in all conditions, identified the most significant eSNP in condition A and then determined whether this eSNP was significant in conditions B,C and D. 19 of 20 such shared eSNPs are tested in GTEx. 6 shared eSNP-eGene pairs are eQTLs in heart tissue (31.6%), and 6 are an eQTL in at least one other tissue (31.6%). In line with our previous eGene analysis, shared eSNPs are more likely to be eSNPs in heart tissue than dynamic eSNPs (Chi-squared test; p=0.01). To further investigate eQTL concordance, we compared the eQTL effect size of our dynamic and shared eQTLs to the eQTL effect size in heart tissue determined by the GTEx consortium. Shared eQTLs overlapping heart tissue eQTLs (n=6) have a higher concordance of effect than dynamic eQTLs overlapping heart tissue eQTLs (n=32; r2=0.73 vs. r2=0.12). We are only able to access the significant eQTLs identified by the GTEx consortium and were therefore unable to determine broader concordance across data sets.

We have replaced the GTEx eGene analysis with this new variant-based analysis in the Results section and Materials and methods section.

2) GWAS overlap – "We next chose to analyze the data using a less stringent gene-based approach using GWAS traits of interest. The second analysis allows for the identification of potentially interesting regions where the causal SNP may be different in the eQTL and GWAS data" I am not sure what this means biologically and what conclusions can be drawn from this. What does it mean if there is no overlap in the SNPs between a eQTL and GWAS loci but that they are both in the same gene region? Is there an enrichment over what would be expected by chance? the statements "We did not mean to imply that the overlap is a shared causal relationship. We have changed the language to reflect this. Our goal is to illustrate that the dynamic eQTLs we identified may play a role in higher-level phenotypes." These sentences seems in opposition to one another.

Our study population are the Yoruba from West Africa. Unfortunately, the vast majority of GWAS studies, including the CAD and MI studies reported here, are carried out in predominantly European populations which have an inherently different underlying LD structure. Thus, while our GWAS analysis includes a SNP based approach, we also felt that casting a wider net connecting eQTLs with genes which have orthogonal evidence of being involved in heart disease was prudent. Specifically, given the underlying differences in LD structure between European and African populations we anticipate that some variants identified by GWAS studies will have weak co-localization with eQTLs identified in our study, but the fact that these SNPs are acting to alter the expression of genes implicated in heart disease suggests they may be important for disease risk.

We have included the caveats of this analysis and explained our approach in the manuscript.

Reviewer #3:

The authors have now addressed most of my questions and concerns. Yet, the result of GWAS colocalisation is still in question. I have previously suggested to perform a proper colocalisation analysis using COLOC package or to report the linkage disequilibrium index (r² value) between GWAS index variant and eQTL lead variant. Unfortunately the authors just provided only P-values for associations which does not support the GWAS locus and eQTL are colocalised. The authors have to expand the cis-window (up to 1Mb) for a handful number of genes overlapping with GWAS loci to perform a proper colocalisation analsysis. It is crucial to provide the posterior probability of colocalisation for such a small number of GWAS overlapping genes.

We appreciate the comments from the reviewer and indeed, failed to include the requested LD analysis in the first revision of this paper. Unfortunately, performing a full colocalization analysis using an extended set of variants that overlap the lead GWAS variant is non-trivial as we cannot easily add new variants to the combined haplotype test without performing a near complete reanalysis of the data. However, we can and have carried out the suggested LD analysis. Using LDlink we tested for the r² value between the lead GWAS variants near EIF2B2 and CARM1 and our eSNPs (in the YRI population). We found the lead GWAS variant associated with CAD risk located near the EIF2B2 gene (rs3832966) was in high LD with our dynamic eSNP (rs12588981) associated with expression of the EIF2B2 gene (r² = 0.96, p < 0.001). Similarly, the lead GWAS variant associated with MI near the CARM1 gene (rs4804142) is in lower, but still significant, LD with our condition B eSNP (rs8105092), which is associated with CARM1 expression (r² 0.32 p < 0.001).

We have added text to describe this analysis in the manuscript in the Results section, and the Materials and methods section. We have also added two supplementary figures of the LD matrix plots (Figure 7—figure supplement 2 and Figure 7—figure supplement 3D).

References:

Giambartolomei, Claudia, Damjan Vukcevic, Eric E. Schadt, Lude Franke, Aroon D. Hingorani, Chris Wallace, and Vincent Plagnol. 2014. "Bayesian Test for Colocalisation between Pairs of Genetic Association Studies Using Summary Statistics." PLoS Genetics 10 (5): e1004383.

Granja, Jeffrey M., M. Ryan Corces, Sarah E. Pierce, S. Tansu Bagdatli, Hani Choudhry, Howard Y. Chang, and William J. Greenleaf. 2020. "ArchR: An Integrative and Scalable Software Package for Single-Cell Chromatin Accessibility Analysis." bioRxiv. https://doi.org/10.1101/2020.04.28.066498.

Kilpinen, Helena, Angela Goncalves, Andreas Leha, Vackar Afzal, Kaur Alasoo, Sofie Ashford, Sendu Bala, et al. 2017. "Common Genetic Variation Drives Molecular Heterogeneity in Human iPSCs." Nature 546 (7658): 370-75.

Author details

1. Michelle C Ward

   1. Department of Medicine, University of Chicago, Chicago, United States
   2. Department of Biochemistry and Molecular Biology, University of Texas Medical Branch, Galveston, United States

   Present address

   Department of Biochemistry and Molecular Biology, University of Texas Medical Branch, Galveston, United States

   Contribution

   Conceptualization, Formal analysis, Investigation, Visualization, Methodology, Writing - original draft, Project administration

   Contributed equally with

   Nicholas E Banovich

   For correspondence

   miward@utmb.edu

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-1485-320X

2. Nicholas E Banovich

   1. Department of Human Genetics, University of Chicago, Chicago, United States
   2. Integrated Cancer Genomics Division, Translational Genomics Research Institute, Phoenix, United States

   Present address

   Integrated Cancer Genomics Division, Translational Genomics Research Institute, Phoenix, United States

   Contribution

   Conceptualization, Formal analysis, Investigation, Methodology, Writing - review and editing

   Contributed equally with

   Michelle C Ward

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-2604-3247

3. Abhishek Sarkar

   Department of Human Genetics, University of Chicago, Chicago, United States

   Contribution

   Formal analysis, Methodology, Writing - review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-4636-9255

4. Matthew Stephens

   1. Department of Human Genetics, University of Chicago, Chicago, United States
   2. Department of Statistics, University of Chicago, Chicago, United States

   Contribution

   Supervision

   Competing interests

   No competing interests declared

5. Yoav Gilad

   1. Department of Medicine, University of Chicago, Chicago, United States
   2. Department of Human Genetics, University of Chicago, Chicago, United States

   Contribution

   Conceptualization, Supervision, Funding acquisition, Writing - review and editing

   For correspondence

   gilad@uchicago.edu

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0001-8284-8926

• 2,875

  Page views

• 284

  Downloads

• 15

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.