DNA methylation patterns from peripheral blood separate coronary artery disease patients with and without heart failure.

AIMS: Natriuretic peptides are useful for diagnosis and prognostication of heart failure of any cause. Now, research aims to discover novel biomarkers that will more specifically define the heart failure phenotype. DNA methylation plays a critical role in the development of cardiovascular disease with the potential to predict fundamental pathogenic processes. There is a lack of data relating DNA methylation in heart failure that specifically focuses on patients with severe multi-vessel coronary artery disease. To begin to address this, we conducted a pilot study uniquely exploring the utility of powerful whole-genome methyl-binding domain-capture sequencing in a cohort of cardiac surgery patients, matched for the severity of their coronary artery disease, aiming to identify candidate peripheral blood DNA methylation markers of ischaemic cardiomyopathy and heart failure. METHODS AND
RESULTS: We recruited a cohort of 20 male patients presenting for coronary artery bypass graft surgery with phenotypic extremes of heart failure but who otherwise share a similar coronary ischaemic burden, age, sex, and ethnicity. Methylation profiling in patient blood samples was performed using methyl-binding domain-capture sequencing. Differentially methylated regions were validated using targeted bisulfite sequencing. Gene set enrichment analysis was performed to identify differences in methylation at or near gene promoters in certain known Reactome pathways. We detected 567 188 methylation peaks of which our general linear model identified 68 significantly differentially methylated regions in heart failure with a false discovery rate <0.05. Of these regions, 48 occurred within gene bodies and 25 were located near enhancer elements, some within coding genes and some in non-coding genes. Gene set enrichment analyses identified 103 significantly enriched gene sets (false discovery rate <0.05) in heart failure. Validation analysis of regions with the strongest differential methylation data was performed for two genes: HDAC9 and the uncharacterized miRNA gene MIR3675. Genes of particular interest as novel candidate markers of the heart failure phenotype with reduced methylation were HDAC9, JARID2, and GREM1 and with increased methylation PDSS2.
CONCLUSIONS: We demonstrate the utility of methyl-binding domain-capture sequencing to evaluate peripheral blood DNA methylation markers in a cohort of cardiac surgical patients with severe multi-vessel coronary artery disease and phenotypic extremes of heart failure. The differential methylation status of specific coding genes identified are candidates for larger longitudinal studies. We have further demonstrated the value and feasibility of examining DNA methylation during the perioperative period to highlight biological pathways and processes contributing to complex phenotypes.

Introduction

Coronary artery disease (CAD) and heart failure (HF) are leading global causes of morbidity and mortality. Natriuretic peptides, such as B‐type natriuretic peptide, and the prohormone, N‐terminal pro‐B‐type natriuretic peptide, are valuable biomarkers utilized for HF diagnosis and prognostication. However, a future goal is to apply an array of biomarkers to detect and guide therapy of a broader spectrum of HF, including in its asymptomatic early stages. Research to discover novel genetic and epigenetic biomarkers associated with myocardial ischaemia, myocardial wall tension/stretch, myocardial inflammation, and fibrosis therefore represents a key area of investigation.

Part of the reparative process following myocardial infarction (MI) requires synthesis of new extracellular matrix that, in some patients, may adversely remodel the left ventricle reducing the ejection fraction during systole. Ischaemic cardiomyopathy does not occur in isolation. Ischaemia and infarction are potent inflammatory triggers associated with local and varying degrees of systemic innate immune activation. A phase of chronic extracellular matrix remodelling involving monocytes and macrophages becomes apparent. The roles of inflammatory cellular subsets in this process are of increasing interest as it becomes apparent that HF may represent a state of dysregulation of inflammatory and fibrotic responses to chronic cardiomyocyte stress. These processes are fundamentally susceptible to regulation at the genomic level, and epigenetic factors, such as DNA methylation, which ultimately determine cellular phenotypes, are likely to play a critical role.

DNA methylation is an epigenetic mark that is increasingly recognized to play a critical role in cardiovascular disease. Notably, methylation variation associated with the development of aortic atheroma is detectable in peripheral blood leucocytes prior to the development of vascular lesions. Furthermore, distinct patterns of peripheral blood DNA methylation have been related to the risk of ischaemic heart disease and coronary events. Variation in cardiomyocyte DNA methylation and gene expression has been demonstrated in end‐stage cardiomyopathy and dilated cardiomyopathy (DCM). , Importantly, significant overlap with the cardiomyocyte DNA methylation changes potentially contributing to cardiomyopathy has been related to changes in peripheral blood DNA methylation. Indeed, these methylation markers of early systolic dysfunction were shown to outperform N‐terminal pro‐B‐type natriuretic peptide. These findings provide strong justification for further studies to identify DNA methylation markers of HF that may precede current clinical predictors and highlight cellular pathways that contribute to the pathogenesis of HF.

There is currently a lack of data focusing specifically on patients with severe multi‐vessel CAD, HF, and peripheral blood DNA methylation. Given the wide‐ranging impact of diabetic and atheromatous processes within the body, it is highly likely that the status of the peripheral blood methylome may indicate important elements that determine the myocardial response to hypoxic, hypertensive, and hyperglycaemic stress. We therefore conducted a pilot methylome‐wide study, as a planned substudy of a large clinical trial, enrolling patients with severe multi‐vessel CAD and angina requiring coronary artery bypass graft surgery. The Aspirin and Tranexamic Acid for Coronary Artery Surgery trial was a large multicentre randomized controlled trial investigating the impact of aspirin and tranexamic acid on bleeding and thrombotic complications following surgery. For this pilot study, we specifically aimed to demonstrate the explorative utility of methyl‐binding domain (MBD)‐capture sequencing to profile the differences in preoperative peripheral blood DNA methylation in a cohort of patients who present with phenotypic extremes of HF but who otherwise share a similar coronary ischaemic burden, age, sex, and ethnicity.

Methods

Patient inclusion

Following ethical approval (The Alfred Hospital Ethics Committee, HREC Number 11/05), we identified 10 men (55–83 years) with multi‐vessel CAD and HF and having elective or semi‐elective coronary artery bypass graft surgery (Table ). These patients were generally (80%) type 2 diabetic. Cases were matched with 10 male diabetic controls (57–81 years) without HF (Table ). Any information suggesting a patient had impaired left ventricular (LV) function or had experienced symptoms of HF resulted in exclusion from the control group. Patients with contributing cardiac valvular pathology (severe mitral regurgitation) were also excluded. Because of the potentially important role of the inflammatory response in HF, patients on corticosteroids or other anti‐inflammatory medications (e.g. COX‐II inhibitors and non‐steroidal anti‐inflammatory medications excluding aspirin) were excluded. Semi‐elective denoted patients who presented with progressive symptoms and signs of HF and/or angina (small ischaemic events such as non‐ST‐elevation MI) and required surgery prior to discharge from hospital. Patients whose HF and/or angina was acutely related to critical coronary events (preoperative ST‐elevation MI, percutaneous coronary intervention on that admission, inotropic support, and intra‐aortic balloon pump support) were excluded. Elective patients were admitted on the day of surgery with stable HF and/or angina.

Table 1

Cases (HF) and controls

Variable	HF (range)	Controls (range)
Number	10	10
Age	67 (55–83)	64 (50–81)
Diabetes	80%	100%
Hypertension	90%	80%
Angina	80%	100%
Ejection fraction (%)	27% (15‐35)	62% (57–71)
HF symptoms	100%	0%
NYHA (1–4)	2.5 (3–4)	1.4 (1–2)
Surgical LV grade (1–4)	3.2 (3–4)	1 (1)
No. of distal coronary grafts	3.7 (3–5)	3.8 (2–6)
% elective vs. semi‐elective	60%	60%

HF, heart failure; LV, left ventricle; NYHA, New York Heart Association.

LV grade: surgical assessment of the left ventricular systolic function. Ejection fraction reported by preoperative transthoracic echocardiography.

Methylation profiling by methyl‐binding domain‐capture sequencing

Whole blood was collected into EDTA tubes, and DNA was extracted using a by DNeasy Blood & Tissue Kit (Qiagen Pty Ltd, Hilden, Germany) and quantified by Qubit fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). DNA was fragmented by sonication. The MethylMiner Methylated DNA Enrichment Kit (Thermo Fisher Scientific, Waltham, MA, USA) was used to enrich methylated DNA fragments from 500 ng of input material. DNA was eluted from beads with 2 M NaCl solution and was purified by NucleoSpin columns to concentrate DNA and remove salts. Five nanograms of methyl‐enriched DNA underwent library construction using the NEBNext Ultra II DNA Library Prep (New England Biolabs, Ipswich, MA, USA), and the barcoded libraries underwent equimolar pooling and were sent to the Australian Genome Research Facility (AGRF, Melbourne, Australia) for 100 cycles, single end sequencing on HiSeq 2500 (Illumina, San Diego, CA, USA) using version four reagents.

Methyl‐binding domain‐capture sequencing data analysis

Reference human genome (GRCh38 primary assembly) was downloaded from Ensembl alongside Ensembl v84 gene annotation in GTF format. Methyl‐binding domain‐capture sequencing (MBD‐seq) reads were quality trimmed using FastX Toolkit (FastX Toolkit 0.0.14) with a minimum base quality threshold of 20 and minimum read length of 20 nt. Reads were mapped to the human genome using BWA‐MEM (Version 0.7.10). CpG dinucleotide content was determined with a custom BASH script. Peak calling was performed using MACS v1.4d software comparing each patient dataset to a pool of inputs we generated previously. Peak size, GC, and CpG proportions were determined with a custom BASH script. Reads mapped to peaks were counted using featureCounts (v1.4.2) with mapping quality ≥10 and otherwise default settings. In parallel, reads were counted at promoters, as defined as regions ±1 kb of transcriptional start sites, generated using BEDTools software. If a gene possessed multiple transcriptional start sites, the read counts were aggregated to gene‐wise counts with featureCounts. Both methylation matrices underwent filtering to discard regions/genes with fewer than 10 reads per sample on average. General linear model feature of edgeR (Version 3.14.0) was used to identify regions/genes with altered methylation in CCF as compared with control after correcting for patient age. Patient age was scaled in R prior to inclusion in the model. CpG island coordinates were downloaded from the UCSC table browser (http://genome.ucsc.edu/cgi‐bin/hgTables). CpG shore coordinates were derived using BEDTools by extracting regions flanking (up to 2 kb) CpG islands. Genome compartments intron, exon, intergenic, gene body, and coding sequence were obtained from the Ensembl v84 gene annotation file or derived with BEDTools. Open chromatin regions were defined by DNase hypersensitivity experiments from the ENCODE project (wgEncodeRegDnaseClusteredV3.bed.gz). Distal regulatory element coordinates were obtained from a previous study, and coordinates were converted to GRCh38 using the liftOver tool (https://genome.ucsc.edu/util.html). Repeat element coordinates were obtained from RepeatMasker data in the UCSC table browser. For pathway analyses, genes were ranked from most up‐methylated to most down‐methylated using the product of negative log2 of the adjusted P‐value and the sign of the fold change as determined by edgeR. Annotation of differentially methylated regions (DMRs) was performed using GREAT for Supporting Information, Table S2, after coordinate conversion from GRCh38 to GRCh19/Hg19 using liftOver. GSEAPreranked (v2.2.2) was used with the ‘classic’ algorithm to determine differential pathways based on promoter methylation from MSigDB v5.

Validation of differentially methylated CpGs with targeted bisulfite sequencing

Profile plots of methylation signal across DMRs were generated using a custom R script. The BiSearch tool was used to design bisulfite PCR primers to amplify DMRs of interest. Primer sequences used are as follows: MIR3675_F:GATTGAAAAAAGTTTAATAAG, MIR3675_R:TCCTAAAAAAACRAAAAAAACC, HDAC9_F:AAATGATTTAAGAGTTAATGAG, and HDAC9_R:AATCCCACCATTAAAATCAAA. Bisulfite conversion of 100 ng DNA was performed using EpiTect Bisulfite Kit (Qiagen). HotStar Tag (Qiagen) was used to amplify regions from bisulfite converted DNA using with 38 cycles of PCR with the recommended thermal cycling. Amplicons were quantified using MultiNA bioanalyzer (Shimadzu Scientific Instruments, Kyoto, Kyoto Japan). Amplicons from the same patient were pooled to equimolar amounts. NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB) was used to generate MiSeq compatible barcoded sequencing libraries. No size selection was performed, and library was amplified with 6 cycles of PCR. Barcoded libraries were combined in equimolar amounts to make a single pooled library. Sequencing was conducted on Illumina MiSeq with the MiSeq Reagent Kit v3 (600 cycles). PEAR was used to assemble forward and reverse reads into a single sequence, which was subsequently mapped to the regions of interest using Bismark. After methylation calling with Bismark methylation extractor, differential methylation analysis was performed using methylKit that uses logistic regression to determine significance using replicate information. Patient age was not considered in the determination of P‐values for differential methylation. Boxplots were generated using the ggplot2 package in the R programming language.

Results

Methyl‐binding domain capture sequencing

Methyl‐binding domain capture was performed on isolated sheared genomic DNA, and these were sequenced to produce an average of 73 million reads per sample [standard deviation (SD) = 8.7 million] (Supporting Information, Table ). Over 99% of these reads mapped to the human genome, and 76.1% (SD = 2.2%) were mapped uniquely (mapQ ≥ 10). As expected, the CpG dinucleotide content in MBD‐seq datasets was elevated (6.54%, SD = 0.65%) when compared with background samples (2.45%). This result indicates that the MBD capture enriched for methylated DNA.

Methylation analysis

Peak calling identified 548k (SD = 55k) enriched regions per sample, occupying 562 Mbp (SD = 65 Mbp) of the genome. Individual peak sets were merged into a list of 643k peaks, which occupied 890 Mbp and had a median size of 973 bp. The CpG dinucleotide content was 3.0% compared with the genome‐wide level of 1.9%. The peak set median GC content was 44.1% compared with the genome‐wide level of 39.3%. We next counted uniquely mapped reads aligned to these regions and omitted peaks below the detection threshold (10 reads per sample on average). This resulted in a matrix of 567 188 detected methylation peaks, accounting for on average 64.3% of MBD‐seq reads (SD = 2.2%). In contrast, only 23.3% of input reads were assigned to peak regions. This result demonstrates that the MBD‐seq identified methylated regions in patients with ischaemic cardiomyopathy with severe multi‐vessel coronary disease requiring coronary artery bypass graft surgery.

Differentially methylated positions in ischaemic cardiomyopathy with heart failure

We implemented a general linear model statistical approach to identify DMRs in HF while correcting for variation in age. We identified 68 DMRs (false discovery rate ≤0.05; Figure ) when compared with 63 DMRs without adjustment for age. As shown in Figure , 48 DMRs occurred within gene bodies and 25 were located near enhancer elements (and full list in Supporting Information, Table S2). Of the DMRs that intersect gene bodies with reduced methylation in HF, 19 were protein‐coding genes (ABCB1, ARID5B, CAV2, DIS3, DOPEY2, GREM1, HDAC9, HS6ST3, IQCA1, JARID2, MYO1H, NFIA, OPA1, PACRG, RPS6KA5, RUNDC3B, SOS2, SOX6, and TC2N) and four non‐coding RNA genes. From DMRs with increased methylation in HF, 19 protein‐coding genes were identified (ACSS2, DAAM2, GPA33/DUSP27, EDEM1, FAM212B, FRMD4A, KCNK10, MICALCL, MYO5A, PARP4, PDSS2, PLA2G4E, RALGAPA1, RGS9, SCAMP1, and TRPM6) and 12 non‐coding RNA genes were identified. To identify broad trends of altered methylation sites across different genomic compartments, we intersected these compartments with DMRs passing a looser significance threshold (nominal P ≤ 0.01). When compared with the genome‐wide background, gene bodies, protein‐coding regions, CpG shores, and regions immediately flanking transcriptional start sites were overrepresented, while DMRs were less frequently identified in CpG islands and intergenic regions (Supporting Information, Figure A). We also observed differential methylation of Alu repeat elements, with Alu Y elements strongly overrepresented in DMRs with elevated methylation in HF, while Alu S elements were overrepresented in DMRs with lower methylation in HF (Supporting Information, Figure B).

Figure 1

Differential methylation among coronary artery disease patients with and without heart failure. (A) Methyl‐binding domain‐capture sequencing identified 567, 188 peaks, of which 68 were differentially methylated (false discovery rate <0.05) between normal and heart failure groups after correction for age, as shown in the volcano plot in red points. (B) Heatmap of the top 50 differentially methylated regions (DMRs) selected by low P‐value.

Gene set enrichment analysis

To further understand the potential functional implications of differential methylation, we performed gene set enrichment analysis (GSEA) for methylation differences at or near known gene promoters in particular known Reactome pathways. We revealed 103 significantly (false discovery rate <0.05) enriched gene sets in HF with promoter down‐methylation occurring more than up‐methylation (87 vs. 16). Table demonstrates top‐ranked selected gene sets with high positive and negative enrichment scores (Table ).

Table 2

Top 10 selected significantly enriched gene sets in coronary artery disease patients with heart failure

ES	Reactome ID (R‐HSA)	Reactome description	n	P‐value	FDR q‐value	Key gene pathway area and function
−0.39	72187	mRNA processing	50	<0.0001	<0.0001	Metabolism of RNA: mRNA splicing and transport
−0.35	159236	mRNA transport	68	<0.0001	<0.0001	Metabolism of RNA: mature RNA transport via nuclear pore complex
−0.31	3214841	Methylation of histones	42	<0.0001	0.0005	Chromatin organization: methylation of lysine 5 of histone H3
−0.59	390666	Serotonin receptors	10	0.002	0.01	Signal transduction: 5‐HT1 and 5‐HT5A receptors that bind serotonin
−0.38	264876	Insulin processing	24	<0.0001	0.01	Metabolism of proteins: translocation of insulin secretory granules
0.30	381070	Activates chaperones	50	<0.0001	0.008	Metabolism of proteins: unfolded protein response in the endoplasmic reticulum
0.42	449836	Other interleukins	23	<0.0001	0.011	Immune system: interleukin signalling: IL‐34, IL‐32, and CD4
0.33	2029481	FC gamma receptor activation	38	<0.0001	0.015	Innate immunity: FCGR‐dependent phagocytosis: immunoglobulin light chains
0.39	5576893	Plateau phase	25	<0.0001	0.018	Cardiac conduction: calcium transport
0.33	445355	Smooth muscle contraction	33	0.02	0.025	Muscle contraction: tropomysin 2 and myosin light chains

ES, enrichment score (negative: reduced promoter methylation; positive: increased promoter methylation); FDR: false discovery rate; n, number of genes in the gene set.

Targeted validation of differential methylation using bisulfite sequencing

Validation of the MBD‐seq data was performed using bisulfite sequencing to determine methylation at two DMRs of the same DNA samples. As the average DMR size was larger than what can be routinely amplified by bisulfite PCR, we visualized methylation signals across DMRs to target groups of CpG nucleotides with the strongest differential methylation. The DMRs selected were located within or nearby to the gene HDAC9 and the uncharacterized miRNA gene MIR3675 (Figure ). While not all targets could be successfully amplified in all samples, at least five individuals per group were obtained for these two regions of interest. High‐throughput 300 bp sequencing of both termini gave coverage of 112k‐fold on average (SD = 148k). The HDAC9 DMR portion analysed, within intron 19 (chr7:18 834 681–18 834 888), had high methylation in the control group (91.8%), and the HF group had lower methylation at two of the four CpG sites analysed. Sites c34 and c49 were lower in HF by 0.62% and 0.99%, respectively (Figure ). The MIR3675 DMR portion analysed (chr1:16 547 840–16 548 225) contained 30 CpG sites, 15 of that had intermediate methylation levels (30–70%) and the remainder were highly methylated (>70%). Overall, MIR3675 methylation levels were 3.7% higher in HF patients, but there were 15 CpG sites where HF patients had methylation 5% higher than the control group (Figure ). These results confirm two regions of altered methylation in HF and validate the use of MBD‐seq for the detection of altered methylation in the perioperative samples.

Figure 2

Targeted bisulfite sequencing‐based validation of two heart failure (HF)‐associated differentially methylated regions (DMRs). (A) Methylation landscape of methyl‐binding domain‐capture sequencing signal across two selected DMRs, HDAC9 and MIR3675. Red arrowheads indicate approximate location on the chromosome ideogram. Black bars indicate CpG sites. Green arrows indicate locations of primers for bisulfite sequencing. Below, tracks of transcription factor binding sites from the ENCODE project data are indicated. (B) Bisulfite sequencing of HDAC9 and MIR3675 that reveal single‐nucleotide‐resolution cytosine methylation at these DMRs depicted in boxplots confirms genome‐wide findings in the same samples.

Discussion

This is the first study utilizing the perioperative environment to create a genomics biorepository to investigate the relationship between peripheral blood DNA methylation and a phenotypic extreme identified within a clinical trial. We demonstrate the utility of MBD‐seq to identify differential methylation at multiple sites throughout the genome in ischaemic cardiomyopathy and HF and to further reveal where altered methylation occurs within different genomic contexts and biological pathways.

Cardiac remodelling and HDAC9

The most well‐described protein‐coding gene identified and validated in this study with links to cardiomyopathy and hypertrophic remodelling is the class IIa histone deacetylase enzyme, HDAC9 (HDAC9; Ch7). HDAC9 is well known in cardiovascular medicine for its action to repress myocyte enhancer factor‐2, a transcription factor that reprograms cardiac gene expression causing cardiac hypertrophy. HDAC9 acts to suppress cardiac hypertrophy, , and our data highlights the methylation status of HDAC9 as a potential marker of a propensity for pathological cardiac remodelling in severe CAD. The data also provide further evidence in support of the therapeutic potential of HDAC inhibitors to treat and potentially prevent the progression to HF. While still at a preclinical stage, HDAC inhibition paradoxically has been shown to be antihypertrophic, the net effect of modified histone acetylation of genes that results in antifibrotic, anti‐inflammatory, and anti‐apoptotic actions on an array of cell types in addition to cardiomyocytes.

Inflammation and HDAC9

Numerous clinical studies have demonstrated an association between reduced vascular events, lipid‐lowering therapies (statins/ezetimibe), and improved outcomes when the level of systemic inflammation is reduced. This is further supported by the observation that lowering systemic inflammation with inhibition of the pro‐inflammatory cytokine, interleukin 1β using canakinumab, in patients with stable CAD reduced cardiovascular mortality by 31%. Recent data for HF in mice implicate a change in the role of myocardial regulatory T lymphocytes (Tregs), becoming pro‐inflammatory, exacerbating inflammation and adverse remodelling. This suggests that Tregs may underlie immune activation in HF and raises another potential effect of HDAC9 activity in this process. By means of its interaction with the transcriptional regulator, forkhead family protein FOXP3 in human Tregs, , , HDAC9 activity and inhibition thereof is a target with widespread immunomodulatory effects that may be of benefit in HF. This hypothesis is supported by the observation that in patients with HF, FOXP3 expression decreases and HDAC9 expression increases as B‐type natriuretic peptide levels rise. The reduced methylation of HDAC9 observed in this study may be an epigenetic marker of increased inflammation mediated by altered Treg activity in HF.

Coenzyme Q10 synthesis and PDSS2

Another highlighted cellular process of interest is the coenzyme Q10 (CoQ10) biosynthetic pathway. This pathway is also well known in HF therapy with dietary supplementation being shown to improve cardiac tolerance to hypoxia, reduce HF symptoms, and reduce major adverse cardiovascular events. , In this study, we identified an increased methylation of the decaprenyl diphosphate synthase subunit 2 (PDSS2; Ch 6) gene. This gene is critical to CoQ10 activity and myocardial mitochondrial function, and the methylation status may prove to be a marker of the capacity of the myocardium to recover from hypoxic stress and potentially indicate therapeutic efficacy of CoQ10 supplementation.

Additional candidate novel biomarkers of heart failure

JARID2 (Ch6) demonstrates reduced methylation. Jumonji (Jumonji and AT‐rich interaction domain containing 2) protein is a transcriptional repressor required for normal murine cardiac development and, with reduced methylation of jarid2, has been associated with greater progression of isoproterenol induced HF. The impact of reduced methylation of JARID2 is unclear; however, in humans, myocardial JARID2 protein is decreased by 60% in HF and altered methylation may actually be an early marker of ventricular decompensation.

GREM1 (Ch15) demonstrates reduced methylation. Gremlin‐1 is an antagonist of bone morphogenic proteins, and its expression has been significantly correlated with myocardial fibrosis and impaired LV ejection fraction in nonischaemic HF. Furthermore, Gremlin‐1 levels have been shown to be elevated in patients with CAD. Its action as an inhibitor of macrophage migration inhibitory factor, leading to increased risk of atherosclerotic plaque instability has led to the proposal that gremlin‐1/macrophage migration inhibitory factor ratio may act as a marker of the severity and instability of CAD.

The enrichment analysis revealed highly ranked genes with decreased promoter methylation involved with mRNA processing (THOC1: Rank 30384), mRNA transport (TPR: Rank 30180), and the nuclear pore complex (NUP54: Rank 30085). Others included diverse areas such as histone methylation (KMT2A: Rank 30305), serotonin receptors (HTR1D: Rank 28393), and myosin VA (MYO5A: Rank 29307), an actin‐based motor protein involved in vesicular transport, insulin processing, and cardiomyocyte potassium channel transport. In contrast, enriched gene sets with increased methylation contained highly ranked genes involved in protein metabolism (TPP1: Rank 421) and components of the nuclear lamina previously associated with DCM (LMNA: Rank 1368). One gene (HYOU1: Rank 1943) is thought to have a cytoprotective role in hypoxia and has been highlighted in DCM. Within the immune system, interleukin signalling (IL32: Rank 284; IL34: Rank 668) and immunoglobulin light chains (IGLV1‐47: Rank 152) were highly ranked. IL34 expression has recently been associated with post‐operative atrial fibrillation. Some genes have been directly related to cardiac conduction (CACNB1: Rank 392, voltage‐gated L‐type calcium channels), contractility (TPM2: Rank 191, tropomyosin 2 slow‐type muscle fibre ; MYL7: Rank 400, myosin light chain 7), and cardiomyopathy.

Overlap with previous studies of DNA methylation and heart failure

There have been previous microarray studies of HF including one by Hass and colleagues using the Illumina 27k Infinium array of LV tissues comprising eight controls and nine patients with DCM, which identified three differentially methylated genes LY75, ERBB3, and ADORA2A. In a larger study of LV (n = 19) and right ventricular (n = 9) tissues in DCM using the 450k Infinium array, 1555 differentially methylated genes were identified, including ARID5B and NFIA, which were also identified here. Although these regions do not overlap, the ARID5B DMR identified here is located nearby three differentially methylated probes (ARID5B +1.2 kb; NFIA +79 and +182 kb). Another microarray study using MeDIP‐chip profiling of a small number of LV tissues identified three DMRs, although these were not significant in the present study. The 2017 study by Meder et al. described mRNA (RNA‐seq) and DNA methylation (Infinium 450k array) in both LV tissue and peripheral blood in DCM, highlighting several probes with robust differences that were replicated in independent cohorts. The Infinium microarray platform measures methylation at 450k and more recently 850k CpG sites, which represents 1.6% and 3.0% of the 28.3 million CpG sites in the genome. This means that 98.4% of genomic CpG sites are yet to be probed in HF. Even with the most comprehensive array available, 97% of CpG sites remain unassessed.

In contrast to previous studies, we selected MBD‐seq over microarrays or reduced representation bisulfite sequencing as this approach has more comprehensive genome coverage. With sufficient read depth, MBD‐seq is able to profile virtually all CpG sites within mappable parts of the genome, containing ~24.3 M CpG sites (86%) albeit without the single base resolution of bisulfite sequencing. To demonstrate, only 13.4% of detected peaks (76 834/570 961) and only 10.1% of DMRs (7/69) overlap a probe with the 450k Infinium array. Similarly, only 5.3% of detected peaks (30 562/570 961) and only 1.4% of DMRs (1/69) overlap a CpG site with ≥10‐fold coverage by ENCODE RRBS (HAIB Methyl RRBS GM12878 Rep1). While MBD‐seq lacks the ability of whole‐genome bisulfite sequencing to call absolute levels of methylation at single‐nucleotide resolution, at the level of genes and DMRs, concordance between the methods is above 90%. Moreover, as MBD‐seq requires fewer reads, it is considerably cheaper than whole‐genome bisulfite sequencing, making it applicable for larger studies. Targeted bisulfite sequencing of a small number of DMRs in the same specimens confirmed the ability of MBD‐seq to identify DMRs, as well as providing absolute methylation quantification of these loci at single‐nucleotide resolution.

Limitations of the current study

This is a pilot study and is not powered or specifically designed to discover HF biomarkers. The study size, male sex, and lack of validation in all patients at all genomic sites preclude conclusions of that nature. However, these preliminary observations demonstrate the utility of this method in perioperative genomic research and provide justification for further, candidate‐based, approaches in a larger cohort of patients, with ongoing co‐measurement of established biomarkers. The complexity of HF, coronary atheroma, and diabetic phenotypes combined with variation in medication administration further limits the generalizability of these findings; however, these data, when combined with other validation cohorts, will contribute to a greater understanding of the relationship between ischaemic cardiomyopathy and peripheral markers of systolic dysfunction.

Conclusions

We have demonstrated the utility of MBD‐seq combined with subsequent targeted bisulphite sequencing to evaluate peripheral blood DNA methylation markers in a cohort of cardiac surgical patients with severe multi‐vessel CAD and phenotypic extremes of HF. The differential methylation status of multiple specific coding genes identified, such as HDAC9, is a candidates for larger longitudinal studies. We have further demonstrated the value and feasibility of examining DNA methylation as part of perioperative medical research to highlight cellular pathways and processes contributing to complex phenotypes within a clinical trial.

Ethics approval and consent to participate

All studies were approved by The Alfred Hospital Ethics Committee (HREC Number: 11/05), and all participants provided informed written consent.

Consent for publication

Not applicable.

Availability of data and material

The dataset supporting the conclusion of this article is available in the NCBI GEO repository (https://www.ncbi.nlm.nih.gov/geo/) under Accession Number GSE134766.

Conflict of interest

None declared.

Author contributions

C.R.B. was responsible for the concept, design, and assembly of the patient cohort. C.R.B., P.S.M., S.W., and K.B. were responsible for the biobank creation. A.W.K., R.T., and H.T. were responsible for the tissue processing and DNA preparation. A.K., I.K., H.K., A.E.‐O., and M.Z. were responsible for next‐generation sequencing and validation analyses. M.Z. was responsible for bioinformatics and data repository creation. C.R.B., M.Z., and K.B. drafted the manuscript. All authors contributed to the critical revision of the manuscript for important intellectual content. All authors read and approved the final manuscript.

Funding

This research was funded by a services support grant through the Baker IDI Heart and Diabetes Institute. C.R.B. is supported by the Australian and New Zealand College of Anaesthetists (ANZCA). K.B. is supported by an E. H. Flack Fellowship. Additional infrastructure funding to the Baker Heart and Diabetes Institute and Murdoch Children's Research Institute was provided by the Australian Government National Health and Medical Research Council Independent Research Institute Infrastructure Support Scheme and the Victorian Government's Operational Infrastructure Support Program. A.E.‐O. is a National Health and Medical Research Council (NHMRC) Senior Research Fellow (1154650 and 0526681). The funding bodies had no role in the design, analysis, interpretation, or reporting of results.

Figure S1. Differential methylation enriched in genomic compartments. Regions with differential methylation (nom p‐val<0.01) were intersected with various annotated functional compartments (A), major repeat element classes (B), and Alu subclasses (C).

Click here for additional data file.

Table S1. Pre‐operative history: Control (GRP 1 CB no 1‐10) Case (GRP 2 CB no 11‐20)

Table S2. Pre‐operative cardiac data: Controls (GRP 1) and Case (GRP 2)

Table S3. Cardiac medications: Control (GRP 1) and Case (GRP 2)

Table S4. Blood tests control: (GRP 1) and Cases (GRP 2)

Click here for additional data file.

Table S5. MBD‐seq data metrics

Table S6. Heart failure DMRs after adjusting for age. DMRs (all FDR≤0.05) are sorted by nominal p‐value. Gene annotations were generated using GREAT. Enhancer annotations were obtained from GeneCards.

Click here for additional data file.