The landscape of DNA methylation modification changes during the progression of AMI
We generated the AMI mouse models by the ligation of the proximal left anterior descending coronary artery, and six time points were selected (Sham, AMI 10-min, 1-h, 6-h, 24-h and 72-h, n = 18) to evaluate the molecular changes in the process of AMI (Fig. 1A). At each time point, the infarcted left ventricular tissues from AMI mice were harvested for RNA-seq and MeDIP-seq. The effectiveness of the mouse models in simulating AMI process has been proved by cardiac ultrasound in the published work by our team [14].
To explore the relationship between DNA methylation modifications and AMI progress, we used the Peaks location information obtained from the MeDIP-seq sequencing data to analyze Peaks, genes and gene functional elements (Promoter, 5’ UTR, CDS, 3’ UTR, Intron and TTR). The CpG islands were annotated to detect the distribution of methylation-enriched regions on the genomic elements (Additional file 1: Figure S1). We found that the overall distribution of Peaks on the genome did not change much in different time points, but it should be some of the minor differences that lead to huge changes in the transcriptional level of genes. In order to obtain the methylation sites related to the development of AMI, we used the sham group as a control to perform differential DNA methylation analysis in different time points post-AMI (|log2FC|> 1; P value < 0.05). As shown in Fig. 1B and Additional file 2, we obtained 18814, 18614, 23587, 26018 and 33788 differential methylation positions (DMPs) at 10-min, 1-h, 6-h, 24-h and 72-h AMI, respectively. And there were 7951(42%), 9108(49%), 1212(51%), 12707 (49%), 13631 (52%) hypomethylation sites of DMPs, respectively (Fig. 1B). From the temporal distribution of DMPs, we found that the numbers and distinctions of DMPs were getting larger with the gradual progress of AMI. Figure 1C shows the top 50 DMPs that appeared simultaneously at these 5 time points.
To visualize the overall expression patterns of DMPs in different time points, we performed the hierarchical clustering analysis of all DMPs around different regions, respectively (Additional file 1: Figure S2). The Promoter and Exon regions showed similar clustering results that no significant difference existed in the methylation expression patterns between 10-min and sham groups, while distinct differences appeared at 1-h, 6-h, 24-h, and changed again at 72-h. The Intron and 3'UTR regions showed that 6-h was clustered into a separate category from 1 h and 24-h, suggesting these two methylation regions did not vary monotonically with time. In addition, the distal intergenic has more independent DNA methylation patterns among different time points, while the clustering effect of the 5'UTR region was not obvious.
A large number of studies have demonstrated that DNA methylation in the promoter region affects gene expression [15], so we specially selected the CpG islands in the promoter regions for further analysis, and identified 2638, 2477, 2979, 3704 and 5881 DMPs at the set time points post-AMI, respectively (Fig. 1D).
Genes differentially expressed become evident in the early AMI phase
To determine the transcriptome alterations in the development of AMI, RNA-Seq was used to detect the mRNA expression in the same sample sets. Our published work has performed the principal component analysis (PCA) on the transcriptome profiles and indicated that the 6-h group, 24-h group and 72-h group are clearly separated from other groups, while the very early time points (10-min, 1-h) seem to be close to the sham group [14]. Here, we identified 123 differentially expressed genes (DEGs) in the 10-min AMI group compared with the control group (|log2FC|> 1 and P value < 0.05), of which 89 genes were upregulated, while the other 34 were downregulated (Fig. 2A). At the following four time points, there were 135, 731, 1419, and 2779 DEGs compared with the control group, respectively (|log2FC|> 1 and P value < 0.05 for 1-h, 6-h and 24-h; |log2FC|> 1.5 and P value < 0.05 for 72-h), and the expressions of most DEGs increased (Fig. 2B–E and Additional file 3). Besides, the results indicated that the number of DEGs and the degree of differences in transcriptional expression gradually increased over time after AMI. It is worth noting that the drastic changes of DEGs started to occur from 6-h post-AMI, while previous stages (10-min, 1-h) showed relatively milder changes.
Dynamic changes of pathways at different time points in AMI
Next, the David gene function annotation online tool was used to perform pathway enrichment analysis on the screened DEGs at each time point. In the very early stage of AMI (10-min to 1-h), a small number of functional pathways were enriched (P < 0.05), including MAPK signaling pathway, TNF signaling pathway and PI3K-Akt signaling pathway, etc. (Additional file 1: Figure S3 A, B). At 6-h AMI, the structure of enriched pathways showed significant alterations that not only lots of new pathways such as Jak-STAT and Adipocytokine signaling pathway were added, but the range of DEGs involved in earlier enriched pathways expanded obviously (Fig. 3). It suggested that 6-h AMI was an important transformation period in the development of AMI. The most significant functional pathways contained TNF signaling pathway, MAPK signaling pathway and cytokine-cytokine receptor interaction, which were mainly associated with inflammation and immune and cellular process. From the stage of 24-h to 72-h AMI, metabolic pathways (Amino sugar and nucleotide sugar metabolism, Protein processing in endoplasmic reticulum) and disease pathways (Hypertrophic cardiomyopathy, Dilated cardiomyopathy) began to add to the enriched functional pathways (Additional file 1: Figure S3 C, D).
To define the temporal characteristics of the complete transcriptome dataset, we performed a clustering analysis of 10,765 genes across 6 time points using Mfuzz, which divided all the genes into 6 clusters (Fig. 4A). The genes in Cluster 2 and 5 both were monotonically upregulated over time, but Cluster 5 rose at 6-h earlier than Cluster 2 at 24-h, while the genes in Cluster 3 gradually downregulated overall. Then, we searched for enriched KEGG pathways (P value < 0.05) within each cluster and summarized these pathways via a heatmap (Fig. 4B). We found that the KEGG pathways were specific to one or more clusters. Some of the pathways enriched in Cluster 3 are unique such as "Oxidative phosphorylation" and "Cardiac muscle contraction," suggesting high temporal specificity with downtrend. On the contrary, there were more overlapping pathways in Cluster 2 and 5, such as "Spliceosome," "Endocytosis" and "Fc gamma R − mediated phagocytosis," indicating an upregulated trend over time.
Identification of DEGs regulated by DMPs in the progression of AMI
Hypermethylation at promoters tends to inhibit the expressions of downstream genes, while hypomethylation increases the expressions [15]. Here, we predicted DEGs regulated by DMPs in AMI process by the following criteria: (1) differential expression of mRNA, (2) differential methylation at promoters, (3) the levels of mRNA and methylation were negatively correlated via Pearson Correlation analysis (P value < 0.05 and Pearson Correlation Coefficient < 0). In the five time points (10-min, 1-h, 6-h, 24-h and 72-h) after AMI, 4, 9, 40, 26, and 183 genes were identified, respectively, that conform to the negative regulatory mechanism of DNA methylation (Fig. 5A and Additional file 4). Figure 5B shows the 40 DEGs regulated by DMPs at 6-h when the large-scale phenotypic changes brought about by DNA methylation started to happen after AMI. Furthermore, we found that most of these identified genes (33 of 40 genes) were highly expressed at transcription level with lower promoter methylation modifications, while the remaining small part of genes were downregulated with higher methylation levels.
Function validation of DNA methylation modification on candidate genes
Based on above results, we performed expression and function validation of DNA methylation on candidate genes that might participate in essential biological process in AMI. First, we manually selected 32 genes which were predicted to be regulated by DNA methylation and involved in important KEGG pathways (Additional file 1: Table S1) to verify their mRNA expressions by qRT-PCR. As shown in Fig. 6A–E and Additional file 1: Figure S4, Spi1, Map3k14, Ncf4, etc., were upregulated at 6 h after MI, Ptpn6, Plcg2, Edem1, etc., were upregulated at 24 h, and Cyba, Itgb5, Col6a1, etc., were upregulated at 72 h.
Then, combining the results of qRT-PCR and MeDIP-Seq data, we excluded genes that the significant changes of DNA methylation occurred after changes in gene expression and then verified the DNA methylation status at promoters of the remaining 21 of above 32 genes using the next-generation sequencing-based BSP (Additional file 1: Table S2). The results showed that the trends in promoter methylation of 10 genes were consistent with predicted outcomes (Fig. 6F–J and Additional file 1: Figure S5). The methylation levels at promoter regions of Ncf4, Map3k14, and Spi1 decreased at 6 h after MI, Ptpn6, Csf1r and Itga11 decreased at 24 h, and Ddost, Cyba, Itgb5, Col6a1 decreased at 72 h.
Furthermore, to address whether the transcription levels of these 10 genes were directly affected by DNA methylation in cardiomyocytes, the whole DNA methylation in H9c2 cells were inhibited in vitro. As shown in Fig. 7, 1 μM Decitabine inhibited the expression of methyltransferase 1 (Dnmt1) in H9c2 cells, and the mRNA levels of Ptpn6, Csf1r, Col6a1, Cyba, and Map3k14 were upregulated after demethylation compared with the untreated group. It indicated that the expressions of these 5 candidate genes were regulated by DNA methylation in cardiomyocytes. Ptpn6, Csf1r, Col6a1, Cyba, and Map3k14 have been involved in the pathogenesis of MI [16,17,18,19,20], so the methylation alterations at their promoters are supposed to play an essential role in AMI.