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

One life-threatening outcome of cardiovascular disease is myocardial infarction, where cardiomyocytes are deprived of oxygen. To study inter-individual differences in response to hypoxia, we established an in vitro model of induced pluripotent stem cell-derived cardiomyocytes from 15 individuals. We measured gene expression levels, chromatin accessibility, and methylation levels in four culturing conditions that correspond to normoxia, hypoxia, and short- or long-term re-oxygenation. We characterized thousands of gene regulatory changes as the cells transition between conditions. Using available genotypes, we identified 1,573 genes with a cis expression quantitative locus (eQTL) in at least one condition, as well as 367 dynamic eQTLs, which are classified as eQTLs in at least one, but not in all conditions. A subset of genes with dynamic eQTLs is associated with complex traits and disease. Our data demonstrate how dynamic genetic effects on gene expression, which are likely relevant for disease, can be uncovered under stress.


Introduction
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, 2019), 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 (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 GTEx project has reported that most eQTLs are shared across tissues, suggesting that they probably do not contribute to disease in a tissue-specific manner.
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 (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;Caliskan et al., 2015;Kim-Hellmuth et al., 2017;Manry et al., 2017;Nedelec 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 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.

Results
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 Methods). iPSC-CMs were matured by electrical pulsing and maintenance in cell culture for 30 days. On Day 30, the median cardiomyocyte purity was 81% (40-97% range), determined by flow cytometry as the proportion of cells that were positive for the cardiac-specific marker, TNNT2 (FigS1; STable; See methods).
We studied the response of the iPSC-CMs to hypoxia and re-oxygenation (Fig1A). To do so, we first cultured the iPSC-CMs at oxygen levels that are close to physiological oxygen levels (10% oxygen -Condition A) for seven days. We then subjected the iPSC-CMs to six 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 (Fig1B,STable). In order to determine whether Figure 1: iPSC-derived cardiomyocytes elicit a cellular response to hypoxia. (A) Experimental design of the study. Cardiomyocytes differentiated from iPSCs (iPSC-CMs) from 15 Yoruba individuals were cultured in normoxic conditions (10% oxygen -condition A) and subjected to 6 hours of hypoxia (1% oxygen -condition B) followed by 6 and 24 hours of re-oxygenation (10% oxygen -conditions C and D). Immunocytochemistry of a representative cardiomyocyte culture where green: TNNT2; blue: nuclei. (B) Peri-cellular oxygen levels of each condition. Each point represents one individual undergoing the oxygen stress experiment. (C) Relative levels of BNP, a marker for cardiac stress, released into cell culture media. Asterisk denotes a statistically significant difference in BNP release (*p < 0.05, **p < 0.005). (D) Molecular phenotypes collected from each individual in each condition. 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.
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), chromatin accessibility data (using ATAC-seq; n = 14), and DNA methylation data (using the EPIC arrays; n = 13; Fig1D) 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 STable and mapped and filtered sequencing reads to prevent allelic mapping biases (FigS3; TableS1; See 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 (FigS4). We filtered out genes with low expression levels (See methods) to yield a final set with data from 12,347 expressed genes (See methods). We performed a number of correlation-based analyses using the data from the technical replicates (FigS5), 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 including MYH7 and TNNT2 (FigS6).
We took advantage of the fact that we have replicate experiments from three individuals to correct the data for unwanted variation (See methods; (Risso et al., 2014)). Following this procedure, our samples clustered both by oxygen level and individual (FigS7). 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;. We used Cormotif (Wei et al., 2015) to classify 2,113 genes (17% of all expressed genes) as responding to hypoxia (Fig2,. Response genes are enriched for genes previously identified to respond to hypoxia in a Caucasian population of individuals (Chi-squared test; P < 2.2x10 -16 ; (Ward and Gilad, 2019)), and are highly enriched for gene ontologies in transcription-related processes (See methods, modified Fisher's exact test; P = 1x10 -19 ). The proportion of all expressed genes that are classified as response genes (green), and non-response genes (magenta). (D) QQ plot illustrating an enrichment of associations between genetic variants and gene expression levels in each condition. Numbers represent the number of eGenes in each condition. (E) An example of a shared eQTL, HEPB2. (F) Heatmap illustrating the 367 SNPs that are classified as dynamic eQTLs. Each row represents a SNP that is an eQTL in at least one condition. Color represents the strength of the association p-value. (G) Examples of each of the two dynamic eQTL categories. Top panel: genes that become an eQTL following hypoxia e.g. ZNF845. Bottom panel: genes that are an eQTL at baseline but not following hypoxia e.g. RFC2. (H) The proportion of baseline eQTLs (all those identified in condition A), and dynamic eQTLs that are also response genes.

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 methods; (van de Geijn et al., 2015)), we identified 1,573 genes with eQTLs (eGenes) in at least one condition (q-value < 0.1; A: 613; B: 564; C: 564 and D: 464; Fig2D). We refer to the 613 eGenes identified in condition A as baseline eGenes. We confirmed that genes whose expression levels are not significantly associated with a SNP in any condition show a largely uniform distribution of P-values in each condition (FigS9A), suggesting that, as expected given our study design, the overall power to detect eQTLs is similar across conditions.
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 naïve comparison of eQTLs classified as 'significant' across conditions will result in an over-estimation of the number of dynamic eQTLs. 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 naïve comparison of eQTLs classified as 'significant' across conditions will result in the identification of true conditionspecific eQTLs). Due to incomplete power, this is clearly not the case (FigS9B); 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 (FigS9B).
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 Pvalue greater than 0.15 (37 instances; Fig2F). 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;Fig2G). This set of 367 dynamic eQTLs corresponds to 328 unique dynamic eGenes (See methods). 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.
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 list of eGenes to eGenes identified in tissue from two locations in the primary heart -left ventricle (LV) and atrial appendage (AA) -from hundreds of individuals in the GTEx study (Consortium et al., 2017). We found that the majority of eGenes present in all conditions ('shared', n=61) in our study are also eGenes in heart tissue (60% in LV, and 67% in AA). A smaller proportion of dynamic eGenes were also classified as primary heart tissue eGenes (38% in both LV and AA). Indeed, shared eGenes are more likely to be eGenes in heart tissue than dynamic eGenes (Chi-squared test; P < 0.005; FigS10A-B). Overall, 98% of dynamic eGenes are an eGene in at least one of 14 selected tissues assayed by the GTEx consortium, supporting our analysis indicating context-dependent inter-individual variation in expression (FigS10C).

Dynamic eGenes are enriched for response genes and transcription factors
To determine whether dynamic eGenes coincide with expression changes of the same genes following hypoxia, we integrated the results of our eQTL and differential expression analyses. In line with our previous findings, baseline eGenes, as well as LV and AA eGenes found in primary tissue (GTEx), are depleted for response genes (Chi-squared test; P < 0.02, (Ward and Gilad, 2019)). However, we found a significant enrichment in response genes amongst dynamic eGenes (61 of 328 genes; Chi-squared test; P = 0.03) when compared to baseline eGenes, suggesting that dynamic eQTLs often impact the regulation of genes that respond to hypoxic stress.
Given that thousands of genes are differentially expressed in response to hypoxia, and many of these genes correspond to transcription-related processes, we next investigated the role of transcription factors, which are likely to drive transcriptome differences in our system. We found a significant enrichment of annotated transcription factors amongst the genes responding to oxygen stress compared to non-response genes (327 of 1,639 annotated human TFs, chi-squared test; P < 2x10 -16 ; FigS11). Given that stress affects transcription factor expression, we asked whether dynamic eGenes are also enriched for transcription factors. Indeed, transcription factors are enriched in dynamic eGenes compared to baseline eGenes (35 TFs; P = 0.004), including MITF and PPARA, both of which are TFs that have been previously implicated in hypoxic response (Feige et al., 2011;Narravula and Colgan, 2001).

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 to identify regions of open chromatin (we were only able to collect these data from 14 of the 15 individuals; STable). We filtered the ATAC-seq reads to include only those reads that map to the nuclear genome, and do not show allelic mapping biases (See methods, 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; however, we identified 831 differentially accessible regions (DARs) between hypoxia and short-term reoxygenation (BC-DARs; 429 regions with increased accessibility and 402 with decreased accessibility), and 71 DARs between hypoxia and long-term re-oxygenation (BD-DARs; Fig3A).
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), and 59 of the 71 BD-DARs are amongst the 831 BC-DARs, suggesting that most regions have returned to baseline levels of accessibility by the first re-oxygenation condition. We therefore considered the 831 BC-DARs, henceforth DARs, in further analysis. Manual inspection of accessibility levels at individual DARs, identified between hypoxia and re-oxygenation, shows that many of these regions appear to have small changes in accessibility between the baseline and hypoxic conditions, which is opposite to the direction of the effect between hypoxia and re-oxygenation. This includes a region within the intron of the FOXO1 gene, a master regulator of the oxidative stress response (Fig3B). Global analysis reveals that there is a strong anti-correlation in the effect size between these pairs of conditions across regions (Spearman correlation = -0.62; sign test P = 4.6 x10 -14 ; FigS15A-C).

Figure 3: Chromatin accessibility changes following hypoxia and re-oxygenation. (A)
Numbers of chromatin regions that are differentially accessible (DARs) between pairs of conditions. (B) Chromatin accessibility levels at a chromatin region within a FOXO1 intron. (C) Expression levels of the hypoxia-responsive gene ADM following hypoxia. (D) Chromatin accessibility levels at a DAR, overlapping an induced HIF1α-bound region, close to the ADM gene.

Linking chromatin accessibility changes with gene expression changes
Changes in chromatin accessibility likely cannot directly explain the thousands of gene expression changes that occur following hypoxia. However, we found that when considering a 50 kb window around the TSS of expressed genes, DARs are enriched near response genes compared to nonresponse genes (Chi-squared test; P = 0.03). 113 of 2,113 response genes have a DAR within 50 kb of the TSS. This set includes an accessible region, overlapping a HIF1α site, within 500 bp of the 3' end of the classic hypoxia response gene, ADM (Fig3C-D).
We asked whether the changes in chromatin accessibility coincide with the appearance of dynamic eQTLs. We found that DARs are no more likely to be near dynamic eGenes than shared or baseline eGenes. In line with previous estimates of the proportion of eQTLs in open chromatin regions, 24 baseline eQTL SNPs (613 total SNPs) and 19 dynamic eQTL SNPs (367 total SNPs) overlap with accessible chromatin regions (Consortium et al., 2017). One dynamic eQTL SNP overlaps a DAR, near the actin filament binding protein gene, FGD4. This gene was also shown to be differentially expressed between children with congenital heart defects where the defect leads to a chronic hypoxic state (cyanotic disease), and children with a similar defect but where oxygen levels are not affected (acyanotic disease; (Ghorbel et al., 2010)).
To directly test whether there are genetic effects on chromatin accessibility, independent of gene expression, we sought to identify chromatin accessibility QTLs (caQTLs) i.e. genetic variants located within the 128,672 accessible regions, which coincide with different levels of accessibility based on genotype. We identified few caQTLs per condition (q-value < 0.1; A: 10, B: 1, C: 7, D: 6; FigS16A). Six of these caQTLs are classified as dynamic caQTLs i.e. induced or suppressed in response to hypoxia using the same definitions as used for the dynamic eQTLs, and include regions at the TSS of the mRNA decapping enzyme gene DCPS, and a region within 100 kb of the C1Orf99 gene (FigS16B-C). These results suggest that gene expression changes, which respond to stress in a genotype-dependent or independent manner, occur largely in the absence of chromatin accessibility changes.

Genomic features associated with differentially accessible regions (DARs)
We next wanted to determine what distinguishes DARs from constitutively accessible regions. 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; P < 2x10 -16 ; Fig4A) 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.2x10 -16 ; H3K4me1: 86% overlap DARs vs. 52% overlap constitutively accessible regions; P < 2.2x 10 -16 ; Fig4B).
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α (Schodel et al., 2011). Based on our inference, 234 of the 356 annotated HIF1α sites (66%), and 150 of the 301 HIF2α sites (50%) overlap with all accessible chromatin regions. We found that HIF1α-and HIF2α-bound regions are more likely to overlap the 831 differentially accessible regions than the 109,275 constitutively accessible regions (Fisher test; P = 0.03; Fig4C).
We next took an unbiased approach to identify transcription factor binding motifs that are enriched in DARs compared to all accessible regions (See methods). We found two motifs to be enriched in DARs compared to all regions (Fig4D). Motif 1 (P = 2x10 -2 ) is recognized by HTF4 and TFE2, both of which are non-response genes in our system. Motif 2 (P = 6.2x10 -42 ) 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 reoxygenation, suggesting that it may be involved in the response. To test this hypothesis, we 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  (Consortium, 2012). (C) The proportion of DARs and CARs that overlap with HIF1α and HIF2α binding locations determined by ChIP-seq in a breast cancer cell line (Schodel et al., 2011). (D) The most significant motif identified to be differentially enriched in DARs compared to all ARs that is putatively recognized by ZNF770, E2F3 and E2F4. We classify E2F4 as a response gene and therefore determined the proportion of DARs and CARs that overlap with E2F4 ChIP-seq binding locations identified in a human LCL line . (E) The proportion of DARs and CARs that overlap four major transposable element (TE) classes -LTR, LINE, SINE, DNA. The proportion of DARs and CARs that overlap with CpG islands (CGIs) that is proximal to the TSS (+/-2 kb from the TSS), and distal to the TSS. Asterisk denotes a statistically significant difference between DARs and CARs (*p < 0.05, **p < 0.005, ***p<0.0005). 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 x10 -6 ; Fig4E).
There is an enrichment of both Alu and MIR SINE family members in DARs (P = 3x10 -5 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.2x10 -16 ; Fig4F).

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 (missing data from two of the individuals), together with 24 replicate samples from three individuals (See methods, STable). 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; FigS17A). Additional analyses indicated the data to be of good quality (FigS17-18).
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 DARs and constitutively accessible regions. To do so, we assessed the DNA methylation level at CpGs within 200 bp upstream of the TSS in the baseline condition. The majority of the assayed CpGs were hypomethylated with a median Beta-value of less than 0.2 across genes and regions (Fig5A).
While there is no difference in median DNA methylation levels between response and nonresponse genes, we found that the median DNA methylation level is lower in DARs compared to all islands. Although we found no difference in CpG methylation between baseline and hypoxia, and the short-term re-oxygenation conditions, we identified four DMCpGs between the baseline and long-term re-oxygenation conditions. This set includes a CpG in the intron of the EGR2 response gene, which shows increased DNA methylation levels over time . Methylation at CpG islands within the intron of EGR2 has been shown to confer enhancer activity in cancer cells (Unoki and Nakamura, 2003). If we only select CpGs located within the promoters of the 2,113 response genes, we find one DMCpG within the promoter of the FTSJ2 gene, a rRNA methyltransferase, that is differentially methylated between hypoxia and long-term re-oxygenation. Selecting CpGs located only within the 831 DARs reveals two DMCpGs between baseline and hypoxia, and one DMCpG between baseline and long-term re-oxygenation. Changes in DNA methylation are therefore not likely to be a major mechanism behind gene expression or chromatin accessibility changes following six hours of hypoxia.

Dynamic eQTLs associate with traits and disease
Finally, we wanted to determine whether any of the dynamic eQTL SNPs or genes that we identified are also associated with complex traits or disease. We first searched within a catalog of genetic variants associated with traits assayed in GWAS for overlap with our dynamic eQTL SNPs (Buniello et al., 2019) and found an induced dynamic eQTL SNPs that is also associated with a measured phenotype -varicose veins (Table1A; Fig6A). We next took an orthogonal approach, using the same GWAS catalog, to specifically investigate three phenotypes that are associated with cardiovascular function or response to oxygen deprivation: MI, heart failure, and stroke (See methods). Six of our dynamic eGenes are also implicated in these disease states by GWAS (Table1B). This list includes the DNA damage and apoptosis factor ZC3HC1, which is implicated in MI and stroke (Fig6B). Importantly, ZC3HC1 is not an eGene in LV or AA, but the SNP-gene pair is an eQTL in other tissues. 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. (A) Dynamic eQTL (deQTL) SNPs associated with GWAS traits. The eQTL SNP, associated gene, type of eQTL based on whether the eQTL is induced or suppressed in response to hypoxia, tested GWAS trait, whether the eQTL SNP is an eSNP in GTEx left ventricle or atrial appendage, and whether it is an eGene in GTEx ventricle or atrial appendage tissue. (B) Dynamic eGenes implicated in three relevant GWAS traits -heart failure, myocardial infarction (MI), and stroke. N.T: not tested.

Discussion
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.

Steady-state and dynamic eQTLs may help understand CVD
Attempts have been made to identify genetic variants that associate with gene expression levels and CVD phenotypes in easily accessible biological samples such as blood. However, less than half of CVD/MI GWAS loci are associated with an eQTL in whole blood when thousands of individuals are tested (Joehanes et al., 2017). To determine the effects of genetic variation on gene expression specifically in the heart, more targeted studies have taken advantage of left ventricle tissue (Consortium et al., 2017;Koopmann et al., 2014), left atrium tissue (Lin et al., 2014;Sigurdsson et al., 2017), and right atrial appendage tissue (Consortium et al., 2017)  there are many loci that remain unexplained.
The heart is a complex tissue consisting of multiple cell types. The effects of some genetic variants in specific cell types might well be masked when considering heterogeneous tissue samples. As we are now able to direct iPSCs towards a cardiac fate, we can test for genetic effects on specific cell types such as cardiomyocytes (Panopoulos et al., 2017). As one would expect, iPSC-CMs are better suited to study cardiovascular traits than the immortalized B-cells or iPSCs from which they are derived (Banovich et al., 2018). However, given the high degree of eQTL sharing across diverse tissues (GTEx), identifying eQTLs in the disease-relevant terminal cell type at steady state may not give substantial insight into disease biology. A significant advantage of using iPSC-CMs is that these cells provide a system to interrogate gene expression dynamics.  (Strober et al., 2019). Further evidence for the notion that steady state eQTLs may have limited applicability to disease states comes from our previous work, where we used a comparative evolutionary approach to investigate the response to stress. We showed that genes that have a conserved response to oxygen stress in iPSC-CMs from both humans and chimpanzees, and are therefore likely relevant for disease, are depleted for eQTLs identified in heart tissue (Ward and Gilad, 2019).
In the current study, by subjecting iPSC-CMs from a panel of individuals to perturbation (oxygen deprivation), we were able to identify a dynamic eQTL SNP (rs8053350) that is associated with varicose veins, and the level of RNF166 expression (Fukaya et al., 2018). This SNP falls within an intron of the PIEZO1 gene. Varicose veins are associated with a risk for developing deep vein thrombosis and other vascular diseases (Chang et al., 2018). When we performed an analogous analysis focused on genes previously associated with three relevant traits -MI, heart failure and stroke, we identified a novel heart eGene, ZC3HC1, encoding the NIPA protein, which is implicated in MI, coronary artery disease and ischemic stroke (Consortium, 2011;Nikpay et al., 2015;Schunkert et al., 2011). This dynamic eQTL SNP is also associated with bronchodilator responsiveness in chronic obstructive pulmonary disease (Hardin et al., 2016).

Mechanisms behind response genes and dynamic eQTLs
Changes in gene expression can associate with other molecular-level phenotypes. The response to hypoxia is mediated by the HIF1α transcription factor (Samanta and Semenza, 2017), but given that there are hundreds of HIF1α binding locations and thousands of differentially expressed genes, regulation by this factor alone cannot directly explain all the transcriptional changes. We Indeed, it has been shown that chromatin contacts exist between HIF1α binding sites and hypoxiainducible genes in the normoxic state (Platt et al., 2016). Conversely, it has been suggested that hypoxia results in the induction of HIF1α, and significant changes in histone methylation (Batie et al., 2019). As we did not measure histone marks in our system, these changes may occur in the absence of chromatin accessibility changes, but we also cannot rule out the possibility that the choice of a single timepoint following six hours of hypoxia, or insufficient statistical power in our sample size, contributed to the minimal differences in accessibility that we observed.
Using an approach designed to measure small effect sizes between conditions, we did identify a set of 831 DARs between hypoxia and short-term re-oxygenation that are enriched for marks of active chromatin, CpG islands, and TEs. These regions do not appear to explain many of the gene expression differences we observed. Hypoxia and oxidative damage are likely to also affect the genome in ways that do not directly impact gene expression. Indeed, the distribution of oxidative DNA damage sites varies across the genome following stress such that TEs and active chromatin regions are enriched for DNA damage, while promoters are depleted (Poetsch et al., 2018). We found enrichment for TEs, specifically Alu SINE elements, in DARs. Interestingly, TEs, and DNA transposons in particular, are also enriched in regions that become accessible in macrophages in response to bacterial infection; suggesting sequence-specific effects of TEs in response to different cellular stressors (Bogdan, 2019). Alu elements have previously been found to associate with the response to stress in other contexts. Serum starvation induces binding of TFIIIC, which recruits RNA polymerase III, to Alu elements (Ferrari et al., 2019), and heat shock increases chromatin accessibility around Alu elements (Kim et al., 2001).
There are several studies, which suggest that DNA methylation levels are dynamic and change in response to stressors such as hypoxia. We did not find any notable differences in DNA methylation levels pre-and post-hypoxia and re-oxygenation, which suggests that like chromatin accessibility, DNA methylation levels do not make large contributions to changes in gene expression levels or the appearance of dynamic eQTLs in our system. Many of the DNA methylation changes that have been described in response to hypoxia occur in chronic and intermittent hypoxia, and not acute hypoxia as investigated in our study (Hartley et al., 2013;Robinson et al., 2012;Watson et al., 2014). DNA methylation levels are also altered in response to other stressors such as bacterial infection ; however, the importance of timing is highlighted by the fact that, in this system, gene expression responses precede DNA methylation changes (Pacis et al., 2019). It is also important to note that our study considers baseline oxygen levels to be 10% oxygen, which is closer to physiological oxygen levels (5-13%) than atmospheric oxygen levels (21%; (Brahimi-Horn and Pouyssegur, 2007;Carreau et al., 2011;Jagannathan et al., 2016). Most studies define normoxia as 21% oxygen saturation, and while this likely leads to larger effect size differences in known hypoxia response genes following hypoxia, these comparisons may not give meaningful insight into the in vivo state.
One can speculate about different mechanisms that might lead to the appearance or disappearance of dynamic eQTLs. In the context of the immune response, it has been shown that the same response variants affect both gene expression and chromatin accessibility (Alasoo et al., 2018). This is in line with the general notion that changing cellular environments results in differences in chromatin accessibility at transcription factor binding sites, which leads to gene expression changes. We found that this does not appear to be a major mechanism in our system as there are minimal changes in accessibility following hypoxia. We observed that there is an enrichment of response genes amongst dynamic eQTLs suggesting that the change in environment results in a change in expression levels that is dependent on the associated genotype. We also find enrichment for TFs amongst response genes and dynamic eQTLs, suggesting that dynamic eQTLs can appear through secondary trans effects.

Potential limitations of our model
To understand the effects of genetic variation on human heart tissue, and how this variation might contribute to the MI and I/R injury etiologies of CVD, we carefully perturbed oxygen levels that cardiomyocytes in culture are exposed to. This in vitro approach is by design a model system, and therefore will likely not fully recapitulate the in vivo state. However, we previously found that out of 2,549 genes that respond to hypoxia in iPSC-CMs from humans and chimpanzees, only 16% are differentially expressed between iPSC-CMs and heart tissue (Pavlovic et al., 2018;Ward and Gilad, 2019). This suggests that our in vitro system is applicable to heart tissue. There is still a possibility that the dynamic eQTLs that we identify in our in vitro system are not physiologically relevant.
Our study comprised a small number of individuals (15), far fewer than what is typical for identifying eQTLs. Our work is therefore a first step towards understanding the effects of genetic variation on gene expression in response to stress. Nevertheless, with a small number of individuals we were able to identify a couple of hundred of dynamic eQTLs that are revealed or suppressed under stress, suggesting that this paradigm is worth exploring further in larger cohorts.
Under the simplifying assumption of a single causal variant, we determined that we have ~6% power to detect an effect which explains ~38% of the heritability, and an equal false positive rate to call it a dynamic eQTL (See methods; FigS20). This suggests that the impact of stress on genotype-dependent effects on gene expression will be even greater in studies which have higher power to detect smaller effects of genotype. For perspective, early eQTL studies were similarly powered to our study, using 70 individuals; yet these studies still led to important insights opening an avenue of research focused on assaying the consequences of genetic variation by RNA-seq (Pickrell et al., 2010).
In summary, there have been few studies assessing the effects of genetic variation in response to CVD-relevant perturbations in cardiomyocytes. Here we profiled the response to oxygen deprivation in cardiomyocytes from a panel of genotyped individuals. We find that eQTLs can appear and disappear in response to oxygen deprivation, and that some of these eQTLs have effects on relevant complex traits and disease.

Cardiomyocyte differentiation from iPSCs
We used fifteen individuals from the Yoruba YRI HapMap population. iPSCs were reprogrammed from lymphoblastoid cell lines (Banovich et al., 2018). iPSCs were maintained in a feeder- On Day 25, iPSC-CMs were transferred to a 10% oxygen environment (representative of in vivo levels) in an oxygen-controlled incubator (HERAcell 150i CO 2 incubator, ThermoFisher Scientific).
From Day 27 onwards, iPSC-CMs were pulsed at a voltage of 6.6 V/cm, frequency of 1 Hz, and pulse frequency of 2 ms using an IonOptix C-Dish & C-Pace EP Culture Pacer to further mature the cells and synchronize beating.

Flow cytometry
Purity of the cardiomyocyte cultures was assessed ~Day 30 as previously described (Ward and Gilad, 2019). Briefly, cells were stained with Zombie Violet Fixable Viability Kit (423113, BioLegend), and PE Mouse Anti-Cardiac Troponin T antibody (564767, clone 13-11, BD Biosciences, San Jose, CA, USA), and analyzed on a BD LSRFortessa Cell Analyzer together with negative control samples of iPSCs, and iPSC-CMs that are incubated without the troponin antibody, or without either the troponin antibody or viability stain.

Hypoxia experiment
On Day 31/32 iPSC-CMs were subjected to the hypoxia experiment. At time = 0, condition A samples remained at 10% O 2 (normoxia), while samples for conditions B, C and D were transferred to an incubator set at 1% O 2 (hypoxia). After 6 hours, conditions A and B were harvested while plates C and D were returned to normoxic oxygen conditions. Plate C was harvested 6 hours following the hypoxic treatment, and Plate D was harvested 24 hours following the hypoxic treatment. Oxygen levels, experienced by the cells in culture, were measured in cultures from each experimental batch using an oxygen-sensitive sensor (SP-PSt3-NAU-D5-YOP, PreSens Precision Sensing GmbH, Regensburg, Germany), optical fiber (NWDV29, Coy, Grass Lake, MI, USA), and oxygen meter (Fibox 3 Transmitter NWDV16, Coy).

Cell culture media for ELISA and cytotoxicity assays
Aliquots of cell culture media from each experiment were centrifuged at 10 000 rpm for 10 min at 4°C to remove cellular debris. The supernatant was stored at -80°C until further use.

Cell pellets for RNA-seq & DNA methylation arrays
Cells from each well of a 6-well plate were washed twice with cold PBS on ice before collection by manual scraping in 1.5 ml PBS. Cells were pelleted by centrifugation at 7 000 rpm for 8 min at 4°C, flash-frozen and stored at -80°C.

ATAC-seq
We performed ATAC-seq in 14 of the 15 individuals we had gene expression data for (STable for details). ATAC-seq libraries were prepared using the Illumina Nextera DNA sample kit. Libraries were amplified for 10-16 cycles depending on the amplification rate of each library. Each library was amplified in a PCR reaction containing 10 µl DNA, 10 ul dH 2 O, 15 µl NMP (PCR master mix), 5 µl PPC (PCR primer cocktail), and 5 µl index N5, 5 µl index N7). PCR conditions were set at 72°C for 5 min, 98°C for 30 sec, 98°C for 10 sec, 63°C for 30 sec, 72°C for 1 min, repeat steps 3-5 4x and hold at 4°C. The number of cycles per library was determined using a qPCR side reaction as described in Buenrostro et al. (Buenrostro et al., 2013). Libraries were purified using Agencourt AMPure XP beads (A63880, Beckman Coulter, IN, USA), and bioanalyzed to determine library quality. 12 or 16 samples were pooled together to generate four master mixes. Each master mix was sequenced 50 bp paired-end on the HiSeq4000 according to the manufacturer's instructions.

Lactate dehydrogenase activity assay
Lactate dehydrogenase activity (LDH) was measured in 5

RNA-seq analysis
Reads were aligned to hg19 using subread align (Liao et al., 2013). The mapped reads were then reprocessed to reduce reference bias for downstream analyses using the WASP pipeline (van de Geijn et al., 2015). Briefly, reads overlapping polymorphisms segregating in our population were remapped to the genome using the true read, and a version of the read with the alternative allele.
Only reads that mapped uniquely to the same locations with both possible alleles were kept. The median number of reads across conditions was similar (A: 34,353,716; B: 33,493,298; C: 33,883,532; D: 38,147,083). The number of filtered reads mapping to genes was quantified using featureCounts within subread (Liao et al., 2014). We obtained measurements for 19,081 genes.
Sample 18852A was an outlier when considering read count correlations between pairs of samples, and was therefore removed prior to subsequent analyses.

Differential expression analysis
We selected autosomal genes for downstream analysis (18,226). Log 2 -transformed counts per million were calculated (Robinson, 2010), and genes with a mean log 2 cpm < 0 were excluded. We used the fact that we have replicate data from three individuals to remove unwanted variation in our data. We used the RUVs function in the RUVSeq package in R (Risso et al., 2014) to identify such factors. By manual inspection, our data segregated by individual or condition after correction with four factors. For the differential expression analysis, we excluded sample replicate one to avoid the outlier sample and randomly selected replicate two, instead of replicate three, for individuals with replicate samples. We used the RUV factors as covariates in our differential expression analysis using the TMM-voom-limma pipeline (Law et al., 2014;Robinson et al., 2010;Smyth, 2004). We used fixed effects for each condition (A, B, C, D), the RUVs factors as covariates, and a random effect for individual, which was implemented using duplicateCorrelation.
Genes with a Benjamini and Hochberg FDR < 0.1 are classified as differentially expressed (Benjamini and Hochberg, 1995).

-Gene expression trajectory analysis
To identify response genes, we used the Cormotif package in R (Wei et al., 2015) to jointly model pairs of tests. We used TMM-normalized log 2 cpm values as input and considered the following pairs of tests: A vs B, B vs. C and B vs. D to determine which genes are changing their expression during the course of the experiment. The best fit was determined to correspond to two correlation motifs or clusters using BIC and AIC. We classified genes as response genes if the probability of differential expression between conditions was > 0.5 in all pairs of tests.

eQTL identification
To map eQTLs, we analyzed the same samples considered in the differential expression analysis.
Given the sample size in this study, we utilized the combined haplotype test (CHT) to identify eQTLs (van de Geijn et al., 2015). This test models both allelic imbalance and total read depth at a region to identify QTLs. We require 50 total counts and 10 ten allele-specific counts for each gene,  (Delaneau et al., 2017;Ongen et al., 2016)). We then computed an adjusted p-value for each SNP-gene pair by taking the CDF of the fitted Beta distribution, evaluated at the reported CHT p-value.
To call significant eQTLs, we estimated q-values for the set of adjusted p-values for each phenotype, and took tests with q < 0.1. The number of eGenes in each condition was determined by taking the most significant SNP-gene association in each condition (i.e. the top SNP). We defined dynamic eQTLs as either: 1) significant only in A (q < 0.1 in A and permutation-adjusted p > 0.15 in B and C and D; suppressed eQTL); 2) significant in at least one of B, C, or D (q < 0.1) and not nominally significant in A (adjusted p > 0.15; induced eQTL).

Power analysis
For QTL mapping, we assume a linear model We defined a dynamic eQTL as either significant only in A, or significant (after Bonferroni correction, in this analysis) in one of B, C, or D and not significant in A. To estimate the false positive rate of dynamic eQTL calling, we asked what was the probability of a SNP passing this definition, assuming the standardized effect size was identical in all four conditions. We then computed phenotypic variance explained, power to detect an eQTL, and false positive rate to call a dynamic eQTL for every choice of standardized effect size .
Overlapping response genes and eGenes with existing gene sets -Gene ontology analysis Gene set enrichment analysis was performed on response genes, and a background set of all expressed genes using the DAVID genomic annotation tool (Huang da et al., 2009a, b). GO Terms related to Biological Processes were selected, and those with a Benjamini-Hochberg controlled FDR < 0.05 were designated as significantly enriched. Each of the five significantly enriched processes relates to transcription ("DNA-templated transcription", "DNA-templated regulation of transcription", "DNA-templated negative regulation of transcription", "negative regulation of transcription from RNA polymerase II promoter", "positive regulation of transcription from RNA polymerase II promoter"). The most significantly enriched GO terms related to Molecular Functions include "transcription factor activity, sequence-specific DNA binding", "nucleic acid binding" and "DNA binding".

-Transcription factors
A list of 1,637 annotated human TFs was obtained from (Lambert et al., 2018), and intersected with our gene sets.
-GTEx eQTLs eQTLs in LV and AA, and twelve other randomly selected tissues (adipose, brain cortex, colon, lung, liver, muscle, pancreas, pituitary, skin, spleen, thyroid, whole blood) were downloaded from v7 in the GTEx portal (www.gtexportal.org). eGenes were selected at 5% FDR in each tissue, and intersected with our gene categories.

ATAC-seq analysis
Paired-end sequencing reads were aligned to hg19 using bowtie2 with default settings (Langmead

Identification of accessible chromatin regions
To generate a unified list of regions with accessible chromatin across conditions and samples, we first used MACS2 (Zhang et al., 2008) to identify peaks within each sample independently. Next, we used BEDtools (Quinlan and Hall, 2010) with the multiIntersectBed function to identify overlapping peaks within each condition separately. Within each condition, we retained peaks with support from more than three individuals and used the mergeBed function to create a conditionspecific consensus. We then combined and merged the bed files across the four conditions to make a final consensus file containing all the filtered accessible regions. The number of reads mapping to accessible chromatin regions was quantified using featureCounts within subread (Liao et al., 2014).

Identification of differentially accessible regions (DARs)
The 128,673 open chromatin regions associated with count data were filtered to include only those regions on the autosomes, and those which had mean log 2 cpm values > 0 for each region. First, to identify differentially accessible regions we used the same limma framework described above for the RNA-seq data. To test for differences between conditions, a linear model with a fixed effect for condition was used together with a random effect for individual. We did not identify any significantly differentially accessible regions with a Benjamini and Hochberg FDR < 0.1. To identify regions with small effect size differences between conditions we used an adaptive shrinkage method implemented in the ashr package in R (Stephens, 2017). We used the regression estimates (regression coefficients, posterior standard errors, and posterior degrees of freedom) generated by limma to calculate a posterior mean (shrunken regression coefficients), FDR, and False Sign Rate (FSR, probability that the sign of the effect size is wrong). We considered regions to be differentially accessible at FSR < 0.1. We denote regions that are not differentially accessible as constitutively accessible regions.

Overlap of DARs with genomic features -TSS
Transcription start sites were obtained from the UCSC

-Transcription factor binding locations
We obtained ChIP-seq data for the hypoxia-responsive factors HIF1α and HIF2α assayed in the MCF-7 breast cancer cell line (Schodel et al., 2011), and E2F4 in the GM06990 lymphoblastoid cell line . Co-ordinates of the 356 HIF1α, 301 HIF2α and 16,245 E2F4 bound regions were converted from hg18 to hg19 using the liftOver tool in the Galaxy platform (http://galaxyproject.org/; (Afgan et al., 2018)).

-Motif enrichment analysis in DARs
We obtained sequences for all accessible regions and differentially accessible regions using the Galaxy platform (Afgan et al., 2018). We used the MEME-ChIP tool within The MEME Suite (Bailey et al., 2009;Machanick and Bailey, 2011) in Differential Enrichment mode to identify motifs differentially enriched in DARs compared to all accessible regions.

-TEs
We obtained repeat annotations from the RepeatMasker track (Jurka, 2000;Smith, 2010) from the UCSC

-CpG islands
We obtained CpG island annotations from the UCSC Table Browser, and overlapped these regions with DARs and constitutively accessible regions.

caQTL identification
The caQTLs were identified in the same manner as described for the eQTLs. However, in the caQTL analyses, we limited tested SNPs to those falling within the peak regions from our consensus file, as opposed to testing variants within 25 kb of the region. Dynamic caQTLs were identified as for dynamic eQTLs.

DNA methylation analysis
To allow for accurate quantification of DNA methylation levels we removed probes overlapping SNPs with a minor allele frequency of > 0.1, and only retained probes with a detection p-value of > 0.75 across samples. Beta-values (ratio of methylated probe intensity and overall probe intensity, and bounded between 0-1) were quantile normalized using lumiN, and, when appropriate, converted to M-values (log 2 ratio of intensities of methylated probe versus unmethylated probe) using lumi (Du et al., 2008).
The methylation level of CpGs coincides with the expected distribution based on their annotated genomic location i.e. low levels of DNA methylation in CpG islands, and higher levels in CpG island shores, and CpG island shelves respectively (FigS17B). Correlation analysis across all pairs of samples, including replicate samples, reveals clustering primarily by individual rather than condition (FigS18).
To measure the DNA methylation level at gene set promoters, we selected CpGs 200 bp upstream of the TSS (TSS200 defined on the array). We considered all CpGs when overlapping with DARs.

Identification of differentially methylated CpGs (DMCpGs)
Differentially methylated CpGs were identified using the same limma framework as described for the RNA-seq data. Analysis was run using both Beta-values and M-values.

Integration with GWAS-implicated genes
We intersected the Reference SNP cluster ID of our dynamic QTLs with the 158,654 SNPs in the NHGRI-EBI GWAS Catalog available from the UCSC We also considered the 'mapped genes' results from GWAS from thee relevant traits: myocardial infarction (EFO_0000612, 89 genes), heart failure (EFO_0003144, 164 genes) and stroke (EFO_0000712, 255 genes), downloaded from the NHGRI-EBI GWAS Catalog in August 2019.
Gene lists were intersected with our response eGenes.

Data Access
All RNA-seq, ATAC-seq and DNA methylation data have been deposited in the Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/) under accession number GSE144426.