This study was approved by the Royal Brompton hospital ethics review board / Brompton and Harefield trust ethics committee (REC approval 10/H0724/18) and is abiding by all the standards of the Declaration of Helsinki. Written informed consent was obtained from the donors prior to their inclusion in the study. Twelve non-diseased valves free from calcification (6 aortic and 6 mitral valves; 10 males: 2 females; age range 42–64 years, mean age 52.2 years, SD 9.9682) were used in this study. After applying inclusion/exclusion criteria three of the twelve valves were excluded from the downstream analysis (Table 1, Supplementary Figure 1). The non-diseased valves were obtained from unused valves of healthy donor hearts, who died of non-cardiac diseases (Table 1). History, macroscopic, and microscopic evaluation were additionally performed to make sure that the donor hearts chosen are free from cardiovascular and valvular complications. The exclusion criteria of donor hearts were previously described in (14).

Aortic and mitral valves from donor hearts were provided by the Royal Brompton Valve biobank. The ischemic time for fresh tissue harvest was set to not exceed 24 h. The aortic and mitral valve leaflets were excised and separately deendothelialized using collagenase II for 10 min at 37°C as the focus of the study was the interstitial cell population. The leaflets were then washed using PBS, snap frozen and stored in −80°C for DNA extraction. DNA was isolated from deendothelialized tissue using the FitAmp™ Blood and Cultured Cell DNA Extraction Kit (Epigentek, NY, USA, catalog #: P-1018) and was subsequently eluted in TE buffer in a total volume of 40 μl. DNA was finally quantified and quality controlled via Qubit fluorescence.

For each sample, 300 ng of DNA was digested for 2 h with the MSP1 enzyme (20U/sample at 37°C) followed by 2 h with TaqαI (20U/sample at 65°C). Digested, CGI enriched DNA fragments <300 bps in length, were selected for using MQ Binding Beads and subsequently collected for bisulfite treatment. Bisulfite treatment was performed using the Methylamp DNA Bisulfite Conversion Kit (Epigentek, NY, USA, catalog #: P-1001). Bisulfite conversion efficiency of the bisulfite-treated DNA was determined by RT-PCR using two pairs of primers against bisulfite-converted DNA (b-actin) and against unconverted DNA (GAPDH), for the same bisulfite-treated DNA samples. Conversion was deemed successful, if more than 99% of the DNA were converted (Epigentek, NY, USA).

Post-bisulfite PBAT-mediated library preparation was performed using the P-1056A kit. First, DNA end polishing and adaptor ligation was performed. This was followed by library amplification using indexed primers and library purification. The final purified library was eluted in 12 μl of water. Assessment of library quality was done via bioanalyzer and KAPA library quantification (Roche, CA, USA). Finally, 10 nM of sample libraries were subjected to single-end enhanced RRBS on Illumina HiSeq 2500.

A summary of the bioinformatics analysis workflow can be found in (Figure 1, Supplementary Figure 2). First, raw reads were subjected to Quality Control (QC) using FastQC version 0.10.1 (15). Trim Galore version 0.3.7 was then used to remove low quality reads, adapters as well as RRBS-related residues that are artificially added during the end-repair step (16). Trimmed reads were then mapped to the UCSC Homo sapiens genome sequence (version hg19) using Bismark version 0.13.0 (17). To permit only up to one mismatch per seed region, the option “-n 1” was set for Bowtie version 1.0.0 utilized by Bismark (18). Methylation information was extracted from Bismark's sorted and filtered mapping results at base resolution using Bismark's methylation extraction software (17). The subsequent analysis was performed in the CpG context. The R package methylKit version 0.9.2 was used for further analysis of the Bismark methylation extraction reports (19). Samples that generated a minimum of 60 million reads (Supplementary Figure 1) were included in the analysis (Table 1) and processed methylKit. Methylation information found in the aforementioned extraction reports was summarized by methylKit over RefSeq promoters, defined as regions located 1 kb before or after a transcription start site (20). Coverage for each promoter was calculated as the sum of the methylated and unmethylated cytosines. The methylation percentage was calculated as a weighted average of cytosine methylation status (21), which corresponds to the overall proportion of methylated cytosines to the sum of all cytosines within the bin (here promoter with the specified coordinates), and which is covered by a minimum number of reads (21). The promoters were subsequently filtered based on coverage (minimal 5 and maximal 99.9 percentile) and merged for comparative analysis with only those promoters that are covered in all replicates being considered. Additional QC steps can be found in the Supplementary Materials Online (Supplementary Figures 3, 4). RRBS data was deposited in EMBL-EBI's European Genome-phenome Archive (EGA) and is accessible through EGA's accession number EGAD00001006303.

The difference in the methylation level of a promoter between aortic and mitral tissue was calculated as the difference between the weighted average of the percent methylation values of all mitral tissue samples and the weighted average of the percent methylation values of all aortic tissue samples at this promoter. A positive methylation difference (hypermethylation) indicated increased methylation in the promoter associated with mitral compared to aortic tissue. The significance of this difference was evaluated by methylKit using logistic regression (19). To correct for multiple hypothesis testing, the sliding linear model (SLIM) method was used by methylKit (19). DM promoters were finally filtered with the cut-off chosen as a methylation difference that is larger than 10 % and a q-value that is <0.05. Annotation of DMRs with genic features, primarily promoters, was carried out by methylKit using the genomation package (19).

The Protein Analysis Through Evolutionary Relationships (PANTHER) tool was utilized to perform functional enrichment analysis with the background of 16,101 promoters. PANTHER was also used to categorize genes whose promoters were found to be DM according to their molecular functions, biological processes, protein classes, pathways and cellular components (22).

Genes associated with DM promoters were used as seed genes to construct PPI networks using NetworkAnalyst (23). The International Molecular Exchange (IMEx) Interactome database, which contains literature-curated comprehensive data from InnateDB, was used by NetworkAnalyst for the generation of the generic network (23). For network construction, methylation difference values of alternative promoters of the same seed gene were replaced by their average before proceeding to network construction (23). The resulting network was trimmed to the minimum to encompass only nodes that connect the original seed genes. Functional enrichment analysis of all seed genes was performed using NetworkAnalyst's Function Explorer using GO, KEGG and Reactome databases. Enrichment p-values were computed based on the hypergeometric test utilized by the Function Explorer (23). NetworkAnalyst's Module Explorer, was used to identify smaller significantly densely connected subnetworks using the Walktrap algorithm (23). Visualization of PANTHER, KEGG and Reactome pathway terms was performed using the GOplot R package (24).

DMR analysis identified 584 significantly DM promoters, of which 305 showed increased methylation in mitral and 279 in aortic valve tissue (Figure 2A, Supplementary Table 1). The gene Repulsive Guidance Molecule A (RGMA) was associated with the most significantly DM promoter (Figures 2B,C, Supplementary Table 1) followed by TBC1 Domain Family Member 32 (TBC1D32) (Figures 2D,E, Supplementary Table 1), B-Cell Lymphoma 3 (BCL3) (Figures 2F,G, Supplementary Table 1) and the long non-coding RNA RP11-1149O23.3 among others (Figure 2A, Supplementary Table 1).

When we investigated genes linked to the significant differentially methylated promoters in isolation, no significant enrichment was found. When functional classification analysis was performed, the genes associated with the DM promoters were grouped into key valve-related functions and pathways (Figures 3A–E, Supplementary Figure 5) with their methylation direction specified in (Supplementary Table 1). Key valve-related pathways such as Wingless/Integrated (WNT)-, Cadherin-, Transforming Growth Factor Beta (TGFB)-, Integrin-, Endothelin-, Platelet-Derived Growth Factor (PDGF)-, NOTCH signaling as well as angiogenesis and general transcription regulation were identified utilizing the PANTHER Pathway database (Figure 3A, Supplementary Table 2).

When the genes were classified by their molecular function, 69.5% of the genes had a binding and catalytic activity including key genes such as Nitric Oxide Synthase 1 (NOS1) (q-value = 0.04, meth.diff = 17.9 %). Other genes exhibited transcription regulator-, translation regulator- and molecular function regulator activities. The latter category included relevant genes such as Apolipoprotein A5 (APOA5) (q-value = 0.037, meth.diff = 15.3%) and Natriuretic Peptide B (NPPB) (q-value = 0.0027, meth.diff = 15.58%) (Figure 3B, Supplementary Table 3).

Among the relevant biological processes categorizing the genes were “biological regulation” and “cellular component organization or biogenesis.” The former included notable genes such as Metallothionein 1F (MT1F) (q-value = 0.004, meth.diff = 17.07%) and the latter included two pertinent ones namely Hyaluronan Synthase 1 (HAS1) (q-value = 0.00015, meth.diff = 22.15%) as well as Actin Alpha Cardiac Muscle 1 (ACTC1) (q-value= 0.049, meth.diff = 11.23%) (Figure 3C, Supplementary Table 4).

Genes were also grouped into relevant protein classes such as gene-specific transcriptional regulator proteins, which included SMAD Family Member 3 (SMAD3) (q-value = 0.0019, meth.diff = 10.25%), cytoskeletal proteins, which contained Myosin Light Chain 2 (MYL2) (q-value = 0.013, meth.diff = 11.14%) and chromatin-binding protein, which encompassed RB Transcriptional Corepressor Like 1 (RBL1) (q-value = 9.5e-06, meth.diff = 22.98%) (Figure 3D, Supplementary Table 5).

Finally, the cellular component categories into which the genes could be classified were membrane- and extracellular regions among others. Relevant genes such as C-Type Lectin Domain Containing 11A (CLEC11A) (q-value = 0.002, meth.diff = 24.72%) were included in the latter and Solute Carrier Family 16 Member 3 (SLC16A3) (q-value = 1.78e-07, meth.diff = 13.4%) in the former category (Figure 3E, Supplementary Table 6).

After functionally categorizing the genes, whose promoters were significantly DM between the aortic and mitral valves on a genome-wide level, we wanted to assess whether we can predict significantly enriched PPIs among those genes as well as the proteins that are necessary to interconnect them to enable a systems-level analysis.

For that purpose, we mapped those DM genes onto a generic PPI database (see Methods) and constructed a PPI network, which ultimately comprised 715 nodes (genes) of which 308 were seed genes and 2,131 edges (PPIs) (Figure 4A, Supplementary Table 7). Among the network's constituent proteins were several hub nodes, the most connected of which were Ubiquitin C (UBC) (Betweenness centrality = 129,377.98; Degree = 200), SMAD3 (Betweenness centrality = 27,152.04; Degree = 75), Ubiquitin Like 4A (UBL4A) (Betweenness centrality = 13,019.92; Degree = 50), Ribosomal Protein S3 (RPS3) (Betweenness centrality = 9,508.1; Degree = 50), Retinoid X Receptor Alpha (RXRA) (Betweenness centrality = 12,661.72; Degree = 47) and SH3 Domain Containing Kinase Binding Protein 1 (SH3KBP1) (Betweenness centrality = 13,005.65; Degree = 44) (Figure 4A, Supplementary Table 7).

In addition to inspecting the topological properties of the constructed network, we performed functional enrichment analysis of both the hyper- and hypomethylated proteins of the subnetwork as well as them connecting these proteins. This analysis showed significant enrichment of pathways and functions.

Many pathways, which are pertinent to valvular mechanisms, were identified via KEGG pathway enrichment analysis including apoptosis, Nuclear Factor Kappa-light-chain-enhancer of activated B cells (NF-κB) signaling pathway, fluid shear stress and atherosclerosis, Tumor necrosis factor (TNF) signaling pathway, osteoclast differentiation, Interleukin 17 (IL-17) signaling pathway, HIF-1 signaling pathway, regulation of actin cytoskeleton, VEGF signaling pathway and TGFB signaling pathway (FDR < 0.05, Table 2).

To investigate whether relevant significantly enriched pathways are present independent of the choice of the pathway database, we furthermore performed Reactome enrichment analysis, which exclusively highlighted additional key valve-related pathways such as immune response, regulation of lipid metabolism, PDGF, NOTCH1, Fibroblast Growth Factor Receptor (FGFR) and transcription (FDR < 0.05, Supplementary Table 8). Similar to KEGG, it showed enrichment for apoptosis-, TFGB-, interleukin- and hypoxia- related pathways (Table 2, Supplementary Table 8).

After inspecting significant valve-related pathways, we checked whether there are enriched biological processes, molecular functions and cellular components related to the molecular regulation of valve-related processes based on GO databases. Relevant biological processes included transcription, apoptosis, growth factor signaling as well as regulation of cellular component organization (FDR < 0.05, Supplementary Table 9). Molecular functions associated with transcription as well as histone acyltransferase activity such as chromatin-, histone deacetylase-, TF-, SMAD- and NF-κB-binding were among the pertinent significantly enriched molecular functions (FDR < 0.05, Supplementary Table 10). Finally, the most relevant significantly enriched cellular components were Transcription Factor II D (TFIID) complex, spliceosomal complex, chromatin, histone deacetylase complex as well as actin cytoskeleton (FDR < 0.05, Supplementary Table 11).

In addition to examining the network structure as a whole, we performed a module analysis (see Methods), to identify subnetworks that show a significantly increased connection density compared to other modules of the parent network. The first identified subnetwork (p-value 3.04e-13) comprised 45 nodes, of which 43 were seed genes and contained the following hub nodes: UBC, also identified as a hub protein of the parent network (Betweenness centrality = 43; Degree = 939) and Ubiquilin 4 (UBQLN4), not identified previously (Betweenness centrality = 3; Degree = 4) (Figure 4B, Supplementary Tables 12, 13).

The second subnetwork (p-value 0.047) contained 12 nodes, of which 5 were seed genes with Breast And Ovarian Cancer Susceptibility Protein 2 (BRCA2) (Betweenness centrality = 33.5; Degree = 5) as well as APOA5 (Betweenness centrality = 32.5; Degree = 4) being its main hub proteins. Interestingly, all of the module's constituent seed nodes exhibited increased methylation in mitral compared to aortic tissue (Figure 4C, Supplementary Table 14). Relevant enriched pathways in this module included homologous recombination (KEGG/Reactome) and Peroxisome Proliferator-Activated Receptor (PPAR) signaling pathway (KEGG) (FDR < 0.05, Supplementary Table 15). The module's biological processes contained lipid homeostasis and cellular response to external stimulus (FDR < 0.05, Supplementary Table 15). Finally, the most relevant enriched molecular functions encompassed DNA-, enzyme- and TF- binding (FDR < 0.05, Supplementary Table 15).

All the utilized analysis methods of genes and pathways show that DNA methylation plays a crucial role in valve development and disease. Indeed, some of the detected genes and pathways are involved in developmental processes such as cardiogenesis, EMT and protein QC and/or diseases such as myxomatous mitral valve (MMV), BAV and AD (Tables 3, 4).

We further generated a GOChord plot to link the pathways in Table 4 to their constituent genes, which are associated with DM promoters (Figure 5A). The plot revealed several pathways sharing the same genes, such as NF-κB signaling, TNF signaling and osteoclast differentiation, which contain TGFB Activated Kinase 1 Binding Protein 1 (TAB1) and Integrin- and PDGF signaling, which include Collagen Type V Alpha 3 (COL5A3) (Figure 5A).

A GOCircle plot visualizes the significance of each pathway in Table 4 and to show the direction of promoter methylation of their constituent genes (Figure 5B). The most significant pathways in Table 4 were NF-κB, TNF signaling and osteoclast differentiation, which share TAB1 (Figure 5A) and are implicated in aortic stenosis (AS) (Table 4).