Epirubicin Alters DNA Methylation Profiles Related to Cardiotoxicity

Background : Epirubicin (EPI) is an important anticancer drug that is well-known for its cardiotoxic side effect. Studying epigenetic modification such as DNA methylation can help to understand the EPI-related toxic mechanisms in cardiac tissue. In this study, we analyzed the DNA methylation profile in a relevant human cell model and inspected the expression of differentially methylated genes at the transcriptome level to understand how changes in DNA methylation could affect gene expression in relation to EPI-induced cardiotoxicity. Methods : Human cardiac microtissues were exposed to either therapeutic or toxic (IC20) EPI doses during 2 weeks. The DNA and RNA were collected from microtissues in triplicates at 2, 8, 24, 72, 168, 240, and 336 hours of exposure. Methylated DNA immunoprecipitation-sequencing (MeDIP-seq) analysis was used to detect DNA methylation levels in EPI-treated and control samples. The MeDIP-seq data were analyzed and processed using the QSEA package with a recently published workflow. RNA sequencing (RNA-seq) was used to measure global gene expression in the same samples. Results : After processing the MeDIP-seq data, we detected 35, 37, 15 candidate genes which show strong methylated alterations between all EPI-treated, EPI therapeutic and EPI toxic dose-treated samples compared to control, respectively. For several genes, gene expressions changed compatibly reflecting the DNA methylation regulation. Conclusions : The observed DNA methylation modifications provide further insights into the EPI-induced cardiotoxicity. Multiple differentially methylated genes under EPI treatment, such as SMARCA4 , PKN1 , RGS12 , DPP9 , NCOR2 , SDHA , POLR2A, and AGPAT3 , have been implicated in different cardiac dysfunction mechanisms. Together with other differentially methylated genes, these genes can be candidates for further investigations of EPI-related toxic mechanisms. Data Repository : The data has been generated by the HeCaToS project (http://www.ebi.ac.uk/biostudies) under accession numbers S-HECA433 and S-HECA434 for the MeDIP-seq data and S-HECA11 for the RNA-seq data. The R code is available on Github (https://github.com/NhanNguyen000/MeDIP).


Introduction
Epirubicin (EPI) is an important anticancer drug that is widely used in multiple types of cancer treatments even though its utilization leads to a high risk of heart failure [1]. To improve the therapeutic application of EPI, clinicians restricted its dose usage, because a very high dose of EPI (around 900 mg/m 2 ) can cause acute heart failure circumstances. However, long-term observational studies have shown that using EPI also at lower doses can still provoke substantial cardiotoxicity [2,3]. Even though researchers have suggested updated signal transduction models [4], and we studied the side effect of EPI on cardiomyocytes on gene expression [5] and protein [6] levels, deeper insights into EPI toxic mechanisms are still on demand. Since epigenetic modification can influence the functional state of genome regions without changing the DNA sequence, epigenetic signals could provide added value to understanding EPI-related toxic mechanisms.
DNA methylation is an influential epigenetic modification, in which a methyl group is added to the fifth carbon position of the cytosine base [7]. DNA methylation is regulated by both methylation and demethylation processes and plays an essential role in different biological activities based on its location in genomic regions. DNA methylation in intergenic regions can repress the expression of potentially harmful genetic elements. DNA methylation in CpG islands or at the first exon can lead to gene silencing [8], while DNA methylations in other gene regions can also be signals for RNA splicing regulators [9]. Research on EPI in gastric cancer demonstrates that the changes in DNA methylation can be beneficial to understanding the biological mechanism of drug resistance [10]. Thus, DNA methylation can also be a useful tool to reveal drug-induced adverse side effects. Studying DNA methylations can thus be allied with gene and protein expression research in order to produce a multi-layered view on EPI-related toxic mechanisms.
Among a wide range of DNA methylation detecting technologies, methylated DNA immunoprecipitationsequencing (MeDIP-seq) has appeared as a cost-effective approach for genome-wide DNA methylation profiling. This method uses a specific antibody to immunoprecipitate methylated DNA; the obtained fractions are then enriched and read by high-throughput sequencing [11]. The outcome of MeDIP-seq is a high dimensional dataset that requires a defined data analysis pipeline to map the se-quencing reads to the reference genome, calculate methylation levels, and identify differentially methylated regions (DMRs) in the genome. While several tools have been developed to accommodate these analysis steps, QSEA is the most recent one that offers a straightforward procedure to inspect MeDIP-seq data [12]. We recently launched a bioinformatics workflow built on QSEA to analyze the DNA methylation status not only on the DMRs level but also on the gene levels [13]. This workflow could refine candidate genes for further investigation and certainly elevates the application of DNA methylation analysis.
In this study, we intently focused on the DNA methylation profiles under EPI exposure as a follow-up to the previous report that has established the aforementioned bioinformatics workflow [13]. By analyzing the genome-wide methylation status of human cardiac tissues exposed to EPI, we identified candidate genes that had strong DNA methylation alterations related to the EPI-induced cardiotoxicity mechanism. We also examined how changes in DNA methylation of those genes affects their expressions on the transcriptome level. Hence, the outcome of this study could suggest new candidate genes with different levels of regulation for EPI-induced cardiotoxicity research.

Cardiac Microtissue Samples
Human cardiac microtissues consisting of 4000 iPSCderived human cardiomyocytes from a female Caucasian donor and 1000 cardiac fibroblasts from a male Caucasian donor were exposed for 2 weeks to EPI either at therapeutic dose or at toxic dose (IC20) levels calculated based on reverse physiologically based pharmacokinetic (PBPK) modeling [14]. The microtissues were collected in triplicates at 2, 8, 24, 72, 168, 240, and 336 hours. EPI was dissolved in 0.1% DMSO before utilization, thus control samples were also exposed to similar DMSO concentrations over time.

MeDIP-seq Data Analysis
After DNA extraction from microtissues, the methylated DNA fragments were isolated by anti-5methylcytosine antibody and then paired-end sequenced (MeDIP-seq) with 50 bp read length [15]. MeDIP-seq paired-end reads were aligned to human reference genome hg38 using Burrows-Wheeler Alignment tool (BWA) version 0.7.17 [16] and converted to .bam files using Samtools version 1.10 [17]. Thereafter, the aligned MeDIP-seq data were processed following the recent bioinformatics workflow for MeDIP-seq data analysis [13] in R version 4.0.5 [18] using the "QSEA", "annotar", "tidyverse", "BSgenome", and "AnnotationDbi" package [12,[19][20][21][22]. The quality of the MeDIP-seq data was reviewed through the quality control functionalities in the QSEA package. We performed the methylation analysis between all EPI-treated and control samples, as well as between either EPI therapeutic or toxic-treated samples compared to control. All the filtering steps were used with the default settings from the bioinformatics workflow [13]: (i) detect DMRs (p-value < 0.01) with pairwise comparisons between EPI-treated and control samples; (ii) select top 5% genes with the highest number of DMRs across their gene regions; (iii) select top 5% genes with the highest number of DMRs in their promoter region, (iv) select genes that had the absolute log2 fold change (Log2FC) ≥0.5. The average p-value and log2 fold change of each gene were calculated based on the average of the p-value and log2 fold change from all DMRs assigned to that gene. The overlapping genes within different DNA methylation analyses were identified using InteractiVenn tools [23]. Gene Ontology (GO) enrichment analysis on differential methylated gene sets was performed with the PANTHER version 14 using the GO molecular function annotation dataset, no correction after Fisher's Exact test, and default reference for Homo sapiens [24].

RNA Sequencing Data Analysis
The RNA in each sample was isolated and measured using Illumina ribosomal RNA-depleted sequencing method with 100 bp paired-end read. Before sequencing, the sample RNA quality and quantity were examined by the Agilent 420 TapeStation and the Qubit™. After sequencing, the sample sequencing quality was checked by FastQC version 0.11.7 [25] and summarized by MultiQC [26]. The RNA sequencing data were then trimmed by Trimmomatic [27] and mapped to the human genome version GRCh38.p12, Ensembl Archive Release 12 93 [28] using RSEM version 1.3.1 [29] and Bowtie2 version 2.3.4.1 [30]. Only samples that had more than 5 million read counts were included for further analysis. The data quality control, preparation, and processing were tangibly reported in the previous study [5]. Due to the limited amount of microtissue after 336 hours exposed to EPI-related toxic dose, only RNA data were available from samples treated with toxic dose at 2, 8, 24, 72, 168, and 240 hours of exposure. Thereafter, the RNA read counts between EPI-treated and control samples were normalized using the "DEseq2" package [31].

Results
The methylation enrichment efficiency was sufficient in all harvested samples, as indicated in the previous report [13]. In general, DNA methylation profiles differed between EPI-treated and control samples (Fig. 1). This indicates that EPI treatment could significantly alter the DNA methylation in cardiac tissues and MeDIP-seq was able to capture these DNA methylation alterations.

DNA Methylation Analysis of All EPI-treated Samples Compared to Control
The DNA methylation analysis between all EPItreated and control samples unveiled 161,356 unique DMRs corresponding to 19,825 genes. Per gene, the high amount  of DMRs, especially in the promoter region, suggests a strong influence of the EPI treatment on the DNA methylation status of that gene compared to control samples. After filtering, 966 genes were in the top 5% genes that had the highest number of DMR regions across gene regions, while 47 out of these 966 genes had the highest number DMRs in the promoter region (Table 1). After selecting genes that had average absolute log2 fold change ≥0.5, the workflow derived 35 genes with strong methylated alterations from the enormous number of detected DMRs (Table 1, Fig. 2).

DNA Methylation Analysis Between EPI-treated and Control Samples Per Dose Exposure
A similar DNA methylation analysis procedure was employed to analyze the DMRs and corresponding differential methylated genes between EPI-treated and control samples per dose. The workflow again distilled the excessive number of detected DMRs into a shortlist of strong differentially methylated genes ( Table 1, Fig. 3). Intriguingly, EPI therapeutic-treated samples showed a higher number of DMRs regions and gradually a higher number of strong differentially methylated genes compared to EPI toxic-treated samples. The DNA methylation analysis between therapeutic-treated samples compared to controls indicated 37 candidates comprising of 15 hyper-methylated and 22 hypo-methylated genes. This is quite comparable with the outcome of DNA methylation analysis between all EPI-treated samples and controls, which had 35 candidates including 9 hyper-methylated and 26 hypomethylated genes. However, the DNA methylation analysis between EPI toxic-treated samples compared to controls demonstrated 19 candidates, of which one gene, SPG7, was hyper-methylated ( Table 1). The GO enrichment analysis demonstrated that the differential methylated genes are involved in different functional classes (Table 2). While a majority of these genes were concentrated in the catalytic activity (GO:0003824) and binding (GO:0005488) groups, the rest engages in regulator and transporter activities.
Furthermore, we identified the overlapping differentially methylated genes within all foregoing DNA methylation analyses. Eight genes were demonstrated as differentially methylated candidate genes after comparing samples treated with either therapeutic or toxic doses compared to controls. However, in the DNA methylation analysis between all EPI-treated samples and control, one of these genes, ATP11A, was not recognized as differentially methylated candidate gene (Fig. 4, Table 3). Interestingly, while SPG7 was hyper-methylated at the EPI toxic-treated conditions (log2FC_avg = 0.84), it was hypo-methylated at EPI therapeutic-treated condition (log2FC_avg = -0.69) and when using all EPI-treated samples compared to controls (log2FC_avg = -0.68) (Supplementary Tables 1-3). By contrast, while the rest of the candidate genes at the EPI toxic-treated conditions were in hypo-methylated status, some of them were in hyper-methylated status at other conditions. For instance, MAD1L1 was hyper-methylated when comparing EPI therapeutic-treated or all EPI-treated samples to control. On the other hand, NCOR2 was hypermethylated at the EPI therapeutic-treated condition but hypo-methylated when comparing all EPI-treated samples or EPI toxic-treated samples to control ( Supplementary  Tables 1-3). Thus, specific doses and how the MeDIP-seq data were processed had influenced on the outcome of the DNA methylation analysis.

From DNA Methylation to Gene Expression
The MeDIP data and the transcriptome data were obtained from the same microtissues exposed to EPI, thus we were able to evaluate the influence of changing DNA methylation status on the gene expression. The gene expression of eight differential methylated genes that were identified for both EPI therapeutic and toxic-treated conditions is shown in Table 3 and Fig. 5. The methylation status of some genes, such as SUN1, demonstrated a coherent relation with the expression on the transcriptome level. SUN1 was hypo-methylated with log2FC_avg = -0.66, -0.79, and -0.59 for the DNA methylation analysis between all EPItreated, EPI therapeutic-treated, and EPI toxic-treated samples versus control, respectively. The expression of SUN1 in almost all EPI-treated samples was higher than its expression in the corresponding control samples (Fig. 5). Nevertheless, for some genes, the regulation at the DNA methylation level could not entirely be related to the gene expression at the transcriptome level. For example, SPG7 was hyper-methylated at the EPI toxic-treated condition and hypo-methylated at EPI therapeutic-treated condition; however, at the transcriptome level, SPG7 was overexpressed in samples treated with both EPI doses compared to its expression in controls (Fig. 5). Similarly, MAD1L1 and NCOR2 were hyper-methylated during EPI therapeutic treatments and hypo-methylated at EPI toxic treatments but their gene expression did not clearly reflect this (Fig. 5). Thus, the change in methylated status could explain the expression regulation of certain genes, but not for all affected genes in the human genome.
Some candidate genes were differential methylated in either EPI therapeutic or toxic-treated conditions, and also demonstrated distinct expressions on the transcriptome in EPI-treated samples from control. For instance, DPP9 was hyper-methylated at the EPI therapeutic-treated condition (log2FC_avg = 0.61, Supplementary Table 2), and its gene expression in EPI-treated samples was also mostly lower than that in control samples (Fig. 6A). While SMARCA4, HDAC4, PKN1, and RGS12 were hypo-methylated at the EPI therapeutic-treated condition (log2FC_avg = -0.50, -   (Fig. 6A). At the EPI toxic-treat condition, SDHA and POLR2A (log2FC_avg = -0.69 and -1.03 respectively, Supplementary Table 3) were hypo-methylated at the DNA methylation level and consequently showed noticeable up-regulation on the transcriptome level (Fig. 6B). By contrast, although AGPAT3 was hypo-methylated at the EPI toxic-treated condition compared to control (log2FC_avg = -0.84, Supplementary Table 3), AGPAT3 had lower RNA expression levels after 24 hours of EPI exposure compared to corresponding control samples (Fig. 6B).

Discussion
EPI is a popular chemotherapeutic agent with cardiotoxic effects. Although we have investigated the impact of EPI on cellular mechanisms on the transcriptome and protein levels [5,6], there is little information about EPI-induced epigenetic modifications. This study demonstrated the undeniable impact of EPI on the DNA methylation profiles in a human cardiac tissue model. Several genes showed strong DNA methylation alterations under EPI treatment. The gene expression data obtained from the same cardiac tissue thereupon disclosed how several EPIinduced changes in DNA methylation could regulate gene expressions on the transcriptome level. Thus, this study certainly provides new insights into the gene expression and regulation at DNA methylated level related to EPI-induced cardiotoxicity.  Samples compared to controls All EPI-treated samples EPI therapeutic-treated samples EPI toxic-treated samples Overlapping differential methylated genes 7 MAD1L1, PRDM15, NCOR2, SUN1, SPG7, ANKRD11, DENND3 Other differential methylated genes 28 30 By using a recently developed MeDIPseq data analysis workflow, we were able to refine the extensive amount of detected DMRs into a shortlist of strongly differentially methylated genes. We detected 35 genes with DNA methylation alterations in cardiac microtissue in vitro exposed to EPI compared to control. When a similar procedure was employed to analyze DNA methylation profile between specific EPI dose treatment and control sample, the outcome in each analysis step (Table 1) resulted in lists of slightly different candidate genes (Table 3). Although there were still overlapped genes between these DNA methylation analyses (Table 3, Fig. 4), it is clear that various sample grouping approaches can generate different outcomes. Furthermore, it also shows that dose-dependence can lead to different numbers of DMRs and corresponding genes. For example, we had 37 differentially methylated genes at the EPI therapeutic-treated condition while we only had 19 differentially methylated genes at the EPI toxic-treated condition (Table 1). These differential methylated genes are involved in molecular interactions, regulations, and transportations (Table 2).
In particular, the change in methylation status of some genes under EPI treatment (Table 3, Figs. 5,6) is potentially associated with EPI cardiotoxic adverse effects. For instance, SMARCA4 (also known as BRG1) was hypomethylated and provoked up-regulation on the transcriptome level in samples treated with EPI therapeutic dose (Fig. 3A, Fig. 6A). SMARCA4 has been known for its critical role in regulating heart muscle development and disease via myosin heavy chain switch. This gene is generally turned off in cardiomyocytes; however, it is re-activated under stress and its level is correlated with hypertrophic cardiomyopathy severity [32]. Thus, EPI could afflict the expression of SMARCA4 via DNA methylation alterations and stimulate cardiac dysfunctions. Similarly, the gene expression of PKN1 was also up-regulated by hypo-methylation in EPI therapeutic-treated conditions (Fig. 3A). A previous study had revealed that the PKN1 activation can initiate cardiac hypertrophy and fibrosis-associated genes expression, and can be involved in heart failure development [33]. RGS12 was hypo-methylated and up-regulated on the transcriptome level in a part of the samples exposed to EPI  due to different times of exposure (Fig. 3A, Fig. 6A). A rodent study has demonstrated that RGS12 contributes to angiotensin II-induced hypertrophy, and its over-expression has been observed in cardiac hypertrophy and heart failure pathology [34]. Another hypo-methylated gene, HDAC4, is notable for rapid histone methylation regulation in re-sponse to elevated cardiac load [35]. On the other hand, DPP9 was hyper-methylated and led to a lower gene expression in EPI therapeutic-treated samples compared to control (Fig. 3A, Fig. 6A). Other research has revealed that drug-induced DPP9 inhibition in cardiomyocytes can impair the CaMKII-PLB and PKC signaling and cause car-diac dysfunction [36]. All these genes were differentially methylated and potentially related to cardiotoxicity even at EPI therapeutic-treated conditions; this is consistent with the inspection of cohort studies in which cancer survivors who underwent EPI treatment have a higher risk of lateonset cardiac disease [2,3].
Eight genes were consistently differently methylated across EPI-treated samples as well as some genes specifically showed strong methylated alterations at the EPI toxictreated condition (Table 3, Figs. 4,5,6). NCOR2 was hypermethylated in the EPI therapeutic-treated condition but hypo-methylated in the EPI toxic-treated conditions, respectively (Fig. 3). As a nuclear receptor, NCOR2 can regulate the expression of other genes and influence the metabolic oxidative balance in cardiomyocytes [37]. A study has suggested that the differential methylation signature of NCOR2 in CD4+ T cells could be a non-invasive biomarker to identify pulmonary arterial hypertension patients [38]. Furthermore, SDHA, POLR2A, and AGPAT3 genes were hypo-methylated at the EPI toxic-treated condition (Fig. 3B) and also play important roles in cardiac dysfunctions. While POLR2A has been considered a stable heart failure reference gene across rodents and humans [39], SDHA participates in the tricarboxylic acid cycle and mitochondrial respiratory chain. The change of SDHA expression, due to the methylation modification at the DNA level (Supplementary Table 3, Fig. 6B), can impact mitochondrial acetyl-CoA homeostasis and energy metabolic which contribute to heart failure [40]. AGPAT3 is also an enzyme involved in mitochondrial oxidation; thereupon, the change in its expression can consequently affect ATP production [41]. Thus, SDHA and AGPAT3 can be potential drivers in the EPI-induced energy metabolic dysregulation and contribute to heart failure development. Some studies have indicated the influence of anthracycline, i.e., epirubicin and doxorubicin, on immune responses such as interleukin-1 (IL-1) or NLRP3 inflammasome and cytokine release [42,43], as well as supplements which have immune-regulating properties that could alleviate the cytotoxicity [44,45]. However, in this study, we did not observed the methylation changes of IL-1 and NLRP3 after analyzing the MeDIP-seq data. It could be that the change in immune response happens at the protein level; and there is no clear change in the DNA methylation level.

Conclusions
This study demonstrated the change of the DNA methylation profile as well as the changing of gene expression as the consequence of DNA methylation alterations under EPI exposure in in vitro human cardiac microtissues. Since epigenetic modification can influence gene expression, differential DNA methylation alterations could offer a supportive explanation to understanding EPI cardiotoxic mechanisms along with transcriptome and proteome study. A handful of genes that had strong EPI-related DNA methy-lation alterations were named as candidates for further investigation. A part of them, such as SMARCA4, PKN1, RGS12, DPP9, NCOR2, SDHA, POLR2A, and AGPAT3, has disclosed their roles in cardiac dysfunctions as well as potential biomarkers for heart failure in different contexts. This is coherent with the well-known EPI cardiotoxicity adverse effects. Those genes, together with other detected candidate genes can be interesting targets for further investigation in EPI-induced toxic mechanisms.

Author Contributions
NN conducted the research, analyzed the data, and drafted the manuscript; RH and ML helped with using the QSEA package. JK and DJ provided suggestions and revised the manuscript. All authors read and approved the final manuscript.

Ethics Approval and Consent to Participate
Not applicable.