Transcriptome-wide map of N6-methyladenosine (m6A) profiling in coronary artery disease (CAD) with clopidogrel resistance

Abstract

Background

Clopidogrel resistance profoundly increases the risk of major cardiovascular events in coronary artery disease (CAD) patients. Here, we comprehensively analyse global m6A modification alterations in clopidogrel-resistant (CR) and non-CR patients.

Methods

After RNA isolation, the RNA transcriptome expression (lncRNA, circRNA, and mRNA) was analysed via RNA-seq, and m6A peaks were identified by MeRIP-seq. The altered m6A methylation sites on mRNAs, lncRNAs, and circRNAs were identified, and then, GO and KEGG pathway analyses were performed. Through joint analysis with RNA-seq and MeRIP-seq data, differentially expressed mRNAs harbouring differentially methylated sites were identified. The changes in m6A regulator levels and the abundance of differentially methylated sites were measured by RT-PCR. The identification of m6A-modified RNAs was verified by m6A-IP-qPCR.

Results

The expression of 2919 hypermethylated and 2519 hypomethylated mRNAs, 192 hypermethylated and 391 hypomethylated lncRNAs, and 375 hypermethylated and 546 hypomethylated circRNAs was shown to be altered in CR patients. The m6A peaks related to CR indicated lower mark density at the CDS region. Functional enrichment analysis revealed that inflammatory pathways and insulin signalling pathways might be involved in the pathological processes underlying CR. The expression of mRNAs (ST5, KDM6B, GLB1L2, and LSM14B), lncRNAs (MSTRG.13776.1 and ENST00000627981.1), and circRNAs (hsa_circ_0070675_CBC1, hsa-circRNA13011-5_CBC1, and hsa-circRNA6406-3_CBC1) was upregulated in CR patients, while the expression of mRNAs (RPS16 and CREG1), lncRNAs (MSTRG.9215.1), and circRNAs (hsa_circ_0082972_CBC1) was downregulated in CR patients. Moreover, m6A regulators (FTO, YTHDF3, and WTAP) were also differentially expressed. An additional combined analysis of gene expression and m6A peaks revealed that the expression of mRNAs (such as ST5, LYPD2, and RPS16 mRNAs) was significantly altered in the CR patients.

Conclusion

The expression of m6A regulators, the RNA transcriptome, and the m6A landscape was altered in CR patients. These findings reveal epitranscriptomic regulation in CR patients, which might be novel therapeutic targets in future.

Introduction

Clopidogrel, a novel antiplatelet drug used to treat CAD patients after percutaneous coronary intervention (PCI), reduces platelet aggregation by inhibiting the activity of the adenosine diphosphate (ADP) receptor P2Y₁₂. However, approximately 10–30% of treated patients still suffer from ischaemic events [1], which might be due to a reduced degree of platelet inhibition. The inability to effectively prevent platelet aggregation, which results in high residual platelet activity, is considered to be indicative of low response to clopidogrel or clopidogrel resistance. It clearly increases the risk of recurrent cardiovascular and cerebrovascular events in CAD patients after PCI.

Complex factors have been reported to influence the clopidogrel response, including clinical factors (such as diabetes), sequence alterations (SNPs, lncRNAs, miRNAs, and mRNAs), and epigenetic modifications (DNA methylation and RNA methylation). Among these factors, RNA methylation is currently a hot research topic. In most eukaryotes, m6A modification in the 5′ cap region of mRNA plays an important role in the maintenance of mRNA stability, variant splicing, polyadenylation, transportation, and translation, and a m6A mark in the 3′ polyA region contributes to the extracellular transport, initial translation, and structural stability of the mRNA along with polyA-binding proteins [2]. The dynamic modification of m6A methylation is related to life processes such as immune responses and stem cell renewal in the body, which can affect the signalling mediated through downstream molecules, cause gene expression imbalance, and alter cell differentiation, homeostasis, and the stress response, leading to the occurrence of diseases or acquisition of pathological states.

However, few studies have reported the effects of m6A modification on the efficacy of cardiovascular drugs. The previous studies have confirmed that METTL3-mediated m6A modification upregulated the expression of the long-chain intergenic noncoding RNA linc00958, leading to increased adipogenesis and affecting the prognosis of hepatocellular carcinoma chemotherapy [3]. The total saponin content in Panax notoginseng can affect the WTAP/p16 signalling pathway through FTO-mediated m6A modification, thereby inhibiting the proliferation, migration, and intimal hyperplasia of vascular smooth muscle cells [4]. However, whether this dynamic modification can affect the antiplatelet effect of clopidogrel is unknown. In the present study, our group sequenced the whole transcriptome from CAD patients, and then, the m6A methylation abundance of the whole transcriptome was analysed. An integrative analysis of RNAs revealed the molecular mechanism underlying clopidogrel resistance.

Results

Characteristics of the patients

There were 46 patients enrolled in our study, and samples from ten patients (five CR and five NC) were used for the systematic analysis of RNA m6A modification abundance, and samples from the remainder of the patients (18 CR and 18 non-CR) were used for validation. The clinical features of the participants are shown in Table 1. Except for that of the platelet function measure, no significant difference in baseline levels was found between the two groups.

Table 1 Characteristics of the patients included in the profiling analysis

A total of 18 CR and 18 non-CR patients with CAD were analysed for the correlation study. The baseline values of the parameters for these patients relisted in Table 2. In addition to the platelet function index and BMI values, the characteristics of the patients in the validation cohort were closely matched.

Table 2 Characteristics of the patients included in the validated cohort

RNA transcriptome expression profiles in CR and non-CR patients

To further explore the relationship between gene expression and m6A modification, we measured the levels of potential markers of heterogeneity after clopidogrel treatment. We characterized the transcriptome in CR and non-CR patients by RNA-seq with high throughput on the basis of circRNAs, lncRNAs, and mRNAs. With thresholds set to |log2FC|> 1 and a P value < 0.05, we identified 583 lncRNAs (391 downregulated and 192 upregulated) and 922 mRNAs (546 downregulated and 376 upregulated). In the light of the expression profile of the circRNA microarray, long probes and short probes were used, and 5448 dysregulated circRNAs (2529 downregulated and 2919 upregulated) were identified at threshold levels of an |log₂FC|> 1 and a P value < 0.05 (Additional file 1: Stable 1, Additional file 2: Stable 2, and Additional file 3: Stable 3).

Overall characteristics of the m6A modification in CR and non-CR patients

Blood samples from ten patients (five CR and five non-CR) were collected and used for MeRIP-seq analysis. As the Venn diagram in Fig. 1A–C displays, 38,071 m⁶A peaks on mRNAs, 21,454 m⁶A peaks of lncRNAs, and 15,944 m⁶A peaks on circRNAs were discovered between CR and non-CR patients. Among these peaks, 751 nonoverlapping m6A peaks were found on mRNAs (Fig. 1A), 446 nonoverlapping m6A peaks were found on lncRNAs (Fig. 1B), and 78 nonoverlapping m6A peaks were found on circRNAs (Fig. 1C) in CR patients. In addition, 114 nonoverlapping m⁶A peaks were found on mRNAs (Fig. 1A), 152 nonoverlapping m⁶A peaks were found on lncRNAs (Fig. 1B), and 45 nonoverlapping m⁶A peaks were found on circRNAs (Fig. 1C) in the non-CR patients.

Fig. 1

Overview of m6A modification peaks on mRNAs, lncRNAs, and circRNAs in CR and non-CR patients. A–C Venn diagram depicting the common and distinct m6A peaks in mRNAs, lncRNAs, and circRNAs between CR and non-CR patients. D–F The numbers of m6A peaks located on mRNA, lncRNA, and circRNA between CR and non-CR patients

Then, we measured the distribution of m⁶A marks in the whole transcriptome. The vast majority of m6A-modified RNAs were covered by fewer than 5 m⁶A peaks, and a small number were covered by more than 5 peaks (Fig. 1E–F).

We used HOMER motif software to analyse the m⁶A peaks. In contrast with that in the CR groups, the motif sequence in the mRNAs was GGACC in the non-CR groups (Fig. 2A), which was consistent with the previously obtained results. To explore the potential function of m6A, we measured the distribution of m6A peaks across the whole transcriptome. We found that the m6A peaks were mostly enriched at the start of the 3′UTR in both the CR and non-CR patients samples (Fig. 2B). m6A peaks in the CR patient samples showed a lower density in the CDS region than those in the non-CR patient samples. In addition, the distribution of m6A peaks in the CR and non-CR patient samples showed had the same tendency, increasing from the 5′UTR to the start of the 3′UTR and decreasing from the start of the 3′UTR to the end of the transcriptome. We determined that the proportion of m6A peaks was the largest in the 5′UTR and the smallest in the intron in the non-CR patient sample (Fig. 2C and D).

Fig. 2

Overview of m6A peaks located on mRNAs, lncRNAs, and circRNAs in CR and non-CR patients. A Motifs enriched of all the identified m6A peaks in CR and non-CR patients. B Diversity of the m6A peak density in the indicated regions between CR and non-CR patients. C Pie chart revealing the regional distributions of m6A peaks in the RNA transcriptome of CR patients. D Pie chart revealing the regional distributions of m6A peaks in the RNA transcriptome of non-CR patients

Differential m6A-modified RNAs in CR patients samples

Differential m6A-modified mRNAs, lncRNAs, and circRNAs in CR patients samples were analysed and compared to those in the non-CR patients samples. Based on an absolute value of log₂FC greater than 1 and a P value less than 0.05, 2919 hypermethylated and 2519 hypomethylated mRNAs were identified, as shown in volcano plots (Fig. 3A). In addition, 192 hypermethylated and 391 hypomethylated lncRNAs (Fig. 3B) and 375 hypermethylated and 546 hypomethylated circRNAs (Fig. 3C) were identified. According to a heatmap (Fig. 4A–C), we found apparent differences in m6A modification peaks among the three kinds of RNAs between CR and non-CR patient samples in a hierarchical clustering analysis. The ten most hypermethylated and most hypomethylated m⁶A peaks on mRNAs, lncRNAs, and circRNAs of CR patient samples are listed in Tables 3 and 4.

Fig. 3

Abnormally m6A-modified mRNAs, lncRNAs, and circRNAs were identified in CR patient samples compared to those in non-CR patients samples. A mRNAs, B lncRNAs, and C circRNAs

Fig. 4

Hierarchical clustering analysis showing the differences in the RNA transcriptome in CR and non-CR patients based on |log₂FC|> 1 and P value < 0.05 criteria. A circRNAs, B lncRNAs, and C mRNAs. In heatmaps, red indicates hypermethylation, and green indicates hypomethylation

Table 3 The ten most hypermethylated m6A modification peaks in CR and non-CR patients

Table 4 The ten most hypomethylated m6A modification peaks in CR and non-CR patients

To measure the distribution of differentially methylated mRNAs, lncRNAs, and circRNAs in chromosomes, we classified differentially hypermethylated and hypomethylated m6A sites by chromosome. We found that the three chromosomes containing the most hypermethylated and hypomethylated m6A sites in mRNAs, lncRNAs, and circRNAs were chromosomes 19 (213), 7 (166), and 17 (152), which are shown in Fig. 5A–C. For mRNAs, the most hypermethylated m6A sites were located on chromosomes 19(129), 17(76), and 16(75); while the most hypomethylated m6A sites were located on chromosomes 19(13), 17(11), and 12(9), as shown in Fig. 5A. For lncRNAs, the most hypermethylated m6A sites were located on chromosomes 7(80), 17(53), and 6(48); while the most hypomethylated m6A sites were located on chromosomes 15(55), 20(18), and 11(13), as shown in Fig. 5B. For circRNAs, the most hypermethylated m6A sites were located on chromosomes 7(10), 19(9), and 5(7); while the most hypomethylated m6A sites were located on chromosomes 6(5), 4(4), and 5(4), as shown in Fig. 5C.

Fig. 5

Differences in the chromosomal distribution of differentially methylated m6A sites between CR and non-CR patients. A–C Chromosomal distribution of differentially methylated m6A sites in three different kinds of RNAs. D–F Statistical analysis showed fold change differences in hypermethylated and hypomethylated m6A sites in transcript regions of mRNA, lncRNA, and circRNA in CR versus non-CR patients

These hypermethylated and hypomethylated m6A sites in mRNA, lncRNA, and circRNA were then classified by transcript region. As shown in Fig. 5D, the fold change was highest for introns with both hypermethylated and hypomethylated m6A sites in mRNA. For the hypermethylated sites in lncRNAs, the fold change was the highest downstream. For the hypomethylated sites in lncRNAs, the fold change was the highest in the promoter (Fig. 5E). Moreover, the fold change was the highest in the promoter for hypermethylated sites in circRNA, and then, it was the highest in the intron for hypomethylated sites in circRNA (Fig. 5F). From all these results, we concluded that there were distinct differences in m6A marks both in terms of their distribution and in their abundance between CR and non-CR patients, indicating that different m6A modification patterns may lead to clopidogrel resistance.

Functional annotation of altered m6A-modified mRNAs, lncRNAs, and circRNAs between CR and non-CR patients

To characterize the potential function of differential m6A methylation patterns between CR and non-CR patients, mRNAs, lncRNAs, and circRNAs with distinct m6A abundance levels were subjected to Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Gene and Genomes pathway analysis.

GO and KEGG enrichment analyses with circRNA-associated genes that were differentially methylated were performed. For the biological process category, altered m6A-modified genes were significantly interrelated with the regulation of the RNA biosynthetic process (Fig. 6A and D). For the cell composition category, hypomethylated circRNA-associated genes were mainly enriched in platelet alpha granule membrane, cytosol, and nuclear part (Fig. 6A), while hypermethylated circRNA-associated genes were significantly associated with the nucleoplasm, nuclear lumen, and nuclear body (Fig. 6D). For the molecular function category, downregulated circRNA-associated genes were involved in ubiquitin-like protein ligase activity, ubiquitin protein ligase activity, and protein binding (Fig. 6A), while upregulated circRNA-associated genes were mainly enriched in nucleic acid binding, transcription coactivator activity, and transcription cofactor activity (Fig. 6D). KEGG pathway analysis indicated that circRNA-associated genes were enriched in inflammatory pathways (such as the MAPK, Wnt, and Notch signalling pathways) (Fig. 7A and B).

Fig. 6

Gene Ontology enrichment analysis of differentially abundant m6A marks between CR and non-CR patients. A circRNAs in CR patients, B lncRNAs in CR patients, C mRNAs in CR patients, D circRNAs in non-CR patients, E lncRNAs in non-CR patients, and F mRNAs in non-CR patients

Fig. 7

KEGG pathway analysis results showing differentially abundant m6A marks between CR and non-CR patients. A circRNAs in CR patients, B circRNAs in non-CR patients, C lncRNAs in CR patients, D lncRNAs in non-CR patients, E mRNAs in CR patients, and F mRNAs in non-CR patients

In Fig. 6B, for the BP, CC, and MF categories, the hypomethylated lncRNA-associated genes were enriched mainly in lipid metabolic process, glycosphingolipid metabolic process, and negative regulation of protein catabolic process; lysosomal lumen, vacuolar lumen, and primary lysosome; hydrolyzing O-glycosyl compounds, protein phosphatase binding, and acting on glycosyl bonds, respectively. The hypermethylated lncRNA-associated genes, as displayed in Fig. 6E, were enriched mainly in UV protection, cAMP catabolic process, and cyclic nucleotide catabolic process (BP); intracellular part, intracellular, and cytoplasm (CC); enzyme binding, protein domain-specific binding, and cAMP binding (MF). KEGG pathway analysis revealed that differentially methylated lncRNA-associated genes were also significantly related to inflammatory pathways (such as Notch and Toll-like receptor) (Fig. 7C and D).

In Fig. 6C, for the biological process category, mRNAs harboured m6A sites with downregulated methylation were involved in cotranslational protein targeting to the membrane, SRP-dependent cotranslational protein targeting to the membrane, and protein targeting to the ER. However, m6A sites with upregulated methylation were highly involved in nucleic acid metabolic processes, transcription, DNA templating, and RNA metabolic processes (Fig. 6F). For CC terms, downregulated m6A site methylation was significantly enriched in cytosolic ribosome, cytosolic part, and ribosomal subunit (Fig. 6C), while upregulated m6A site methylation was enriched in the nucleus, nuclear lumen, and nucleoplasm (Fig. 6F). Finally, these mRNAs were significantly correlated with nucleic acid binding, an MF category term (Fig. 6C and F). Furthermore, these downregulated mRNAs were found to be significantly involved in the insulin signalling pathway, VEGF signalling pathway, and so on. The upregulated mRNAs were also involved in the insulin signalling pathway and some inflammatory pathways (such as Notch and mTOR) (Fig. 7E and F).

Comprehensive analysis of MeRIP-seq and mRNA-seq data in CR and non-CR patients

The 25 most highly differentially expressed mRNAs are displayed in Table 5. Through a comprehensive analysis of the data from RNA-seq and MeRIP-seq of CR and NCR patients, we identified one hypermethylated and upregulated (hyper-up) gene (LYPD2), one hypomethylated and downregulated (hypo-down) gene (RPS16), and two hypermethylated and downregulated (hyper-down) genes (KDM6B and ST5) with differential expression in CR patients compared to that in non-CR patients (Fig. 8).

Table 5 The 25 most differentially expressed mRNAs based on P value

Fig. 8

Comprehensive analysis of mRNA m6A peaks and their expression between CR and non-CR patients. Scatter plot showing the distribution of mRNAs with both significantly changed m6A and mRNA levels in CR and non-CR patients

Expression of m6A regulators and candidate genes correlated with clopidogrel resistance

In an effort to determine whether m6A methylation status was changed in CR patients, with the validated cohort data, we assessed m6A regulator expression in CR and non-CR patients through RT‒PCR. The results indicated that the expression of FTO, a main m6A eraser (P = 0.0003), and YTHDF3, a main m6A reader (P = 0.0106), was significantly upregulated in CR patients (Fig. 9F and L) compared with non-CR patients (controls). However, the expression of WTAP, a main m6A writer, was lower (P = 0.0244) (Fig. 9D). The expression of other m6A regulators was unchanged between the CR and non-CR patient groups (Fig. 9A–C, E, and G–K).

Fig. 9

The mRNA expression of m6A regulators in CR and non-CR patients. A–E m6A writers (METTL3, METTL14, METTL4, WTAP, and KIAA1429), F, G m6A erasers (FTO and ALKBH5), and H–L m6A readers (YTHDC1-2 as well as YTHDF1-3)

To validate our m6A-seq results, they were performed gene-specific m6A-IP qPCR assays of four differentially expressed genes identified after a comprehensive analysis (KDM6B, ST5, LYPD2, and RPS16), three hypermethylated mRNAs (PIP5K1C, C14orf142, and S1PR2), and three hypomethylated mRNAs (RPS21, SRRT, and RPS6); three hypermethylated lncRNAs (ENST00000584377 (CSNK1D), ENST00000576490 (RP11-473M20.14), and ENST00000438116 (RPL4P6)) and three hypomethylated lncRNAs (ENST00000532093 (THAP1), ENST00000559291 (C15orf57), and ENST00000593731 (PIK3R2)); and three hypermethylated circRNAs (hsa_circ_0096507 (RSF1), hsa_circ_0095501 (NUCB2), and hsa_circ_0138956 (TMEM2)) and three hypomethylated circRNAs (hsa_circ_0007094 (NDRG3), hsa_circ_0104702 (HOMER2), and hsa_circ_0103595 (TP53BP1)) that might participate in pathological processes in CR patients. We observed almost the same relationships between m6A abundance changes and the expression of these genes (except KDM6B, C14orf142, S1PR2, ENST00000438116, ENST00000593731, and hsa_circ_0104702), confirming the validity of our MeRIP-seq results (Fig. 10A).

Fig. 10

Validation of the m6A-enriched genes and candidate loci. A Validation of the m6A-enriched genes as indicated by m6A-immunoprecipitation (IP)-qPCR. B Validation of the expression levels of candidate mRNAs, lncRNAs, and circRNAs as determined by RT-qPCR

Sequentially, to confirm our RNA-seq results, transcript levels of the abovementioned genes (KDM6B, ST5, LYPD2, and RPS16), three upregulated mRNAs (GLB1L2, LSM14B, and CCDC26), lncRNAs (MSTRG.13776.1, ENST00000429224.1, and ENST00000627981.1), and circRNAs (hsa_circ_0070675_CBC1, hsa-circRNA13011-5_CBC1, and hsa-circRNA6406-3_CBC1), as well as downregulated mRNAs (IGLC3, CREG1, and LAMP5), lncRNAs (MSTRG.9215.1, MSTRG.13920.3, and ENST00000508832.2), and circRNAs (hsa_circ_0112493_CBC1, hsa-circRNA393-4_CBC1, and hsa_circ_0082972_CBC1) were also measured in 18 pairs of CR and non-CR patients (Fig. 10B). These findings confirmed that the expression of most of these genes showed similar to the tendency in transcript level changes that we found with the RNA-seq data obtained from the original five pairs of samples.

Discussion

Antiplatelet treatment is the key treatment for patients with CAD after PCI. Although ticagrelor, a new antiplatelet drug, exhibits rapid, effective, and sustained inhibition of the P2Y12 receptor, patient haemorrhage risk is markedly increased [5]. In clinical practice in China, although ticagrelor reduces cardiovascular mortality and the risk of stent thrombosis in ACS patients, it has an advantage only in patients with low bleeding risk and increases the incidence of cardiovascular events in patients with medium to high bleeding risk [6]. The 2018 COSTIC study proposed that clopidogrel shows a significantly greater clinical benefit for the Chinese population [7]. For patients with stable coronary heart disease or elderly patients, clopidogrel might be more beneficial.

However, clopidogrel resistance significantly increases the risk of adverse cardiovascular and cerebrovascular events. The relevant mechanisms underlying the occurrence of clopidogrel resistance are complex. Various studies have indicated that the extent of clopidogrel resistance is influenced by a series of clinical morbidities, including chronic kidney disease, diabetes mellitus, and smoking and factors such as drug–drug interactions (for instance, proton-pump inhibitors), age, inflammatory status, decreased left ventricular ejection fraction, and increased body mass index [8]. The molecular mechanism underlying clopidogrel resistance has been explored in multiple fields on the basis of sequence differences and epigenetic regulatory factors. First, based on single-nucleotide polymorphisms (SNPs) in drug metabolism-related genes, researchers found that ABCB1 genes are involved in duodenal transport; CYP family genes are involved in liver metabolism; PON-1 genes and P2Y12 genes are related to platelet activation. Preliminary research by this project team revealed that the ABCB1 C3435T polymorphism is associated with platelet activity and MACEs in ACS patients taking clopidogrel [9]. Another genome-wide association study suggested that genetic variations in genes such as CYP2C19, PEAR1, and N6AMT1 are associated with different antiplatelet effects and clinical endpoint events related to clopidogrel [10]. Second, the up- or downregulation of noncoding RNA also affects the antiplatelet activity of clopidogrel. A recent omics study suggested that upregulation of the lncRNAs NONHSAT083775.2 and NONHSAT107804.2 or downregulation of the lncRNA NONHSAT133455.2 led to clopidogrel resistance [11]. Another study with DM animal models revealed that lncRNA MT1P3 (metallothionein 1 pseudogene 3) promoted the expression of P2Y12 receptors by sponging miR-126, thereby increasing platelet activity [12]. Moreover, miRNAs participate in the regulation of clopidogrel resistance [13]; for instance, miRNA-26a upregulation increased the clopidogrel resistance risk in patients undergoing PCI [14], and platelet-derived miRNAs (such as miRNA-24-3p, miRNA-142-3p, and miRNA-411-3p) might be potential biomarkers. In addition, the advancement of epigenetic studies has provided a new perspective on the clopidogrel resistance mechanism. Our group found that in CAD patients with hyperlipidaemia, CpG4 methylation of the PON1 promoter resulted in low mRNA expression of PON1, leading to clopidogrel resistance [15]. Our latest research was based on a ChIP analysis of methylated DNA from clopidogrel resistance samples and revealed the identification of 979 hypermethylation sites and 6119 hypomethylation sites, and four sites (cg23371584, cg15971518, cg04481923, and cg22507406) were verified to be clopidogrel resistance diagnostic markers [16]. In addition, N6-methyladenosine (m6A) might play a vital role in posttranscriptional modification.

In the present study, by preparing the RNA transcriptome expression profiles and validating the findings with samples from a larger cohort, we discovered that the expression of mRNAs (ST5, KDM6B, GLB1L2, and LSM14B), lncRNAs (MSTRG.13776.1 and ENST00000627981.1), and circRNAs (hsa_circ_0070675_CBC1, hsa-circRNA13011-5_CBC1, and hsa-circRNA6406-3_CBC1) was upregulated in CR patients, while mRNAs (RPS16 and CREG1), lncRNAs (MSTRG.9215.1), and circRNAs (hsa_circ_0082972_CBC1) were downregulated in CR patients. This is the first study to report on the landscape of the transcriptome based on sequencing of samples from patients with clopidogrel resistance, and it revealed some CR-related biomarkers. Although the results need to be validated with data from a larger sample and from a broader population, the markers it revealed will provide a basis for early detection in clinical practice. Thus, the ST5, KDM6B, and CREG1 genes deserve our continued attention. First, the Cryptococcus neoformans sequence type 5 (ST5) lineage can be used to infect immunocompetent hosts and increase the blood platelet count [17]. In our study, the ST5 level was higher in CR patients, which might have been due to a high level of ST5 increasing the platelet count, leading to clopidogrel resistance. Second, lysine-specific demethylase 6B (KDM6b) is involved in vascular remodelling, and it can promote neointima formation by inhibiting the activation of Nox4 signalling in autophagy [18]. Third, the cellular repressor of E1A-stimulated genes (CREG1), which was downregulated in CR patients, promotes HUVEC proliferation and prevents endothelial cells (ECs) from undergoing apoptosis. Through ILK-Cdc42 activation, CREG1 can also increase EC filopodia formation and vascular assembly to promote neovascularization [19]. Likewise, a latest study uncovered the role of CREG1 in the regulation of megakaryocyte maturation and thrombopoiesis [20]. Hence, ST5, KDM6B, and CREG1 might be therapeutic targets for clopidogrel resistance and vascular diseases. Future researches might further validate the findings of this article.

Moreover, we carried out MeRIP-seq with samples from patients with low clopidogrel responses and found that m6A modification of mRNA on ST5, LYPD2, RPS16, PIP5K1C, RPS21, SRRT, and RPS6 affected mRNA expression, resulting in clopidogrel resistance. In addition, the m6A regulators FTO, YTHDF3, and WTAP were also differentially expressed in CR and non-CR patients. These findings may indicate that upregulated FTO or downregulated WTAP causes reduced expression of related genes. FTO, a key m6A enzyme, is involved in the pathophysiological processes leading to myocardial hypertrophy and cardiac dysfunction [21]. Depending on the m6A modification, FTO regulates the cardiac function after systolic heart dysfunction in model rats and influences cardiac remodelling and repair [22]. Loss of the m6A reader YTHDF2 promotes demethylation of H3K27me3, leading to enhanced production of proinflammatory cytokines [23]. Moreover, WTAP has been reported to regulate the progression of myocardial ischaemia/reperfusion injury via the YTHDF1/FOXO3a signalling pathway [24]. After our comprehensive analysis, we speculated that the higher degree of m6A methylation on ST5 mRNA might lower its expression and result in clopidogrel resistance. The intrinsic molecular mechanism needs to be confirmed with cell experiments.

Furthermore, in this study, we proposed that inflammatory pathways (such as the NF-κB, Wnt, and Notch signalling pathways) as well as the insulin signalling pathway participate in the regulation of clopidogrel resistance. For example, through the NF-κB signalling pathway, nucleotide-binding oligomerization domain 2 (NOD2) mediates P2Y12 upregulation and increases platelet activation, resulting in clopidogrel resistance [25]. On the other hand, abnormal excretion of insulin leads to an increase in platelet activity, and hyperglycaemia activates the protein kinase C pathway by saccharifying proteins on the surface of platelets, inducing the expression of P selectin, reducing the fluidity of the cell membrane, and thus increasing the platelet adhesion rate and enhancing platelet activity [26]. FTO is an important factor in blood glucose regulation, because it regulates insulin signalling-related gene expression possibly through m6A modification. Research suggested that FTO is associated with insulin secretion and abnormal blood glucose levels, and it may interact with transcription factors (STAT3, CREB, ATF4, and FoxO1) and thus change their activity, regulating the expression of liver gluconeogenesis genes (PCK1 and G6PC) [27]. A functional experiment with pancreatic ß cells showed that increasing the expression of FTO led to an increase in the gluconeogenesis rate and that reducing the activity of FTO led to the opposite results [28]. Hence, FTO may regulate insulin signalling-related gene expression through m6A modification, changing the antiplatelet effect of clopidogrel. Further animal and cell experiments are needed to confirm this possibility.

To the best of our knowledge, our research is the first to report RNA transcriptome expression profiles and potential impact of m6A methylation on clopidogrel resistance. Although we made great efforts, this research still has unavoidable limitations. First, it was a single-centre study, and the studied population was not large. We will further expand the sample size to verify our conclusions. Second, PBMCs are collected from peripheral blood, and their status may be influenced by various factors, such as individual differences, etc., which may interfere with the interpretation of RNA transcriptome, and RNA methylation may be influenced by other factors, such as environment factors and coexisting diseases. These influences may have affected our results. Finally, knockout animal models need to be used for in vivo validation studies, and an in-depth mechanistic investigation will increase our understanding in the future.

Conclusions

In general, our present study revealed the RNA transcriptome expression profiles in CR and non-CR patients and showed that the expression of some mRNAs (ST5, KDM6B, GLB1L2, and LSM14B), lncRNAs (MSTRG.13776.1 and ENST00000627981.1), and circRNAs (hsa_circ_0070675_CBC1, hsa-circRNA13011-5_CBC1, and hsa-circRNA6406-3_CBC1) was upregulated in CR patients, while that of other mRNAs (RPS16 and CREG1), lncRNAs (MSTRG.9215.1), and circRNAs (hsa_circ_0082972_CBC1) was downregulated in CR patients. Moreover, m6A regulators (FTO, YTHDF3, and WTAP) were differentially expressed. By MeRIP-seq, we found that 2919 hypermethylated and 2519 hypomethylated mRNAs, 192 hypermethylated and 391 hypomethylated lncRNAs, and 375 hypermethylated and 546 hypomethylated circRNAs were markedly changed in CR patient samples. The joint analysis of mRNA m6A peaks and expression revealed that the mRNAs (such as ST5, LYPD2, and RPS16) might cause clopidogrel resistance because their expression is changed via m6A modification. These are critical regulators that interfere with the epigenetic regulation of clopidogrel resistance via the inflammatory pathway and insulin signalling pathway, and these findings provide new insights useful for CR mechanistic research and will help clinicians identify individualized treatments more effectively.

Materials and methods

Patients

Between April and December 2020, a total of 23 CR and 23 non-CR controls in the First Affiliated Hospital of Ningbo University were enrolled in the present research. All the enrolled patients were diagnosed with ACS and were treated of PCI with drug-eluting stents to re-establish the coronary artery blood supply. The case group conformed to the standard for a CR diagnosis, which is reaction units of P2Y12 (PRU) greater than 240 [29] as measured by the VerifyNow P2Y12 assay (Accumetrics, Inc., San Diego, California). The enrolled CR and non-CR patients did not have a history of acute or severe infection, disordered kidney or hepatic function, active haemorrhage, or rheumatoid-related diseases.

The clinical data, demographic data, and laboratory indices of these patients were collected. These data included BMI, age, sex, drinking rate, smoking rate, comorbidities (diabetes, hypertension, or hyperlipidaemia), medications, aTnI, ALT, AST, LDL-C, creatinine, uric acid, platelets, pro-BNP status, and other measures.

This research was approved by the Ethics Committee of the First Affiliated Hospital of Ningbo University, and every individual signed a written informed consent form before the study was initiated. The study was performed in accordance with Declaration of Helsinki principles.

RNA preparation

Peripheral blood was collected from 23 CR patients and an equal number of non-CR subjects after clopidogrel therapy that was administered longer than 6 months. Samples from 10 of these patients were used for a comprehensive analysis of RNA m6A modification (the five pairs of CR and non-CR patients were well matched), while samples from the remainder of the patients (18 CR and 18 non-CR patients) were used to validate the results.

Peripheral samples were obtained from the patients described above and placed into an ETDA anticoagulant vacutainer. Peripheral blood mononuclear cells (PBMCs) were isolated as described previously [30]. Total RNA was obtained with TRIzol reagent (Invitrogen, USA), and the RNA purity and concentration were determined with a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA).

RT-PCR

RT-PCR was performed with cDNA. First, reverse transcription was performed with RNA by SuperScript III Reverse Transcriptase (Invitrogen, USA) according to the manufacturer’s instructions. The following reagents were added in sequence to a 0.5-ml PCR tube: 2 µl of first-strand cDNA; 2 µl of upstream primer (10 pM); 2 µl of downstream primer (10 pM); 4 µl of dNTP (2 mM); 5 µl of 10 × PCR buffer; and 1 µl of Taq enzyme (2 u/ul). Then, PCR was performed according to the following programme: 95 °C, 10 min; 30 PCR cycles (95 °C, 10 s and 60 °C, 60 s), and amplification reaction 28–32 cycles (95 °C, 10 s; 60 °C, 60 s; and 95 °C, 15 s) and slowly warmed from 60 to 99 °C to establish the melting curve. The housekeeping genes U6/GAPDH were utilized as the internal controls to ensure the reliability and accuracy of the experimental results.

RNA-seq and MeRIP-seq

RNA-seq and MeRIP-seq were performed at CloudSeq Biotech Inc. according to a procedure described previously [31]. Briefly, after RNA isolation, the total RNA was fragmented into approximate 100 nt sequences, and those fragments were incubated with anti-m6A antibody (Manga) for 2 h at 4 °C. Then, total RNA was incubated with prepared beads for 2 h at 4 °C. Finally, the complexes were washed from the beads and purified with TE buffer. Both the RNA that was not immunoprecipitated and IP samples were used for library construction. The library quality was evaluated with a Bioanalyzer 2100 system. The library was sequenced on an Illumina HiSeq instrument. All the sequence data have been uploaded Supplement date (Additional file 4: Stable 4, Additional file 5: Stable 5, and Additional file 6: Stable 6).

Data analysis of MeRIP-seq and RNA-seq

Paired-end reads were obtained using a Illumina HiSeq 4000 sequencer, and the quality was controlled on the basis of a Q30 sore. Then, 3′ adaptor trimming and low-quality reads were removed by Cutadapt software (version: 1.9.3). Then, all the clean reads in the library were aligned to the reference genome (UCSC HG19) by Hisat2 software (version: 2.0.4). DCC software was used to recognize circRNAs based on the results of STAR alignment. Methylated sites on RNAs (m6A peaks) were recognized using MACs software. The recognized methylated sites were used for motif enrichment analysis by HOMER [32]. In addition, the R package MetaPlotR was used to characterize the m6A distribution. diffReps was used to identify differentially methylated sites at threshold levels of |log₂ FC| higher than 1 and a P value less than 0.05. Scripts were written in-house to identify and select the differentially methylated sites that overlapped with mRNA exons. For RNA-seq, raw counts were obtained with HTSeq software (version 0.8.2), followed by normalization using edge R software. Genes of interest were directly visualized with Integrative Genomics Viewer (IGV; version: 2.3.68).

GO and KEGG pathway enrichment analysis

Gene Ontology (GO) enrichment analysis in three categories: biological process (BP), cellular component (CC), and molecular function (MF) categories, was performed to identify the function of identified genes related to differential peaks or differential expression of mRNAs, lncRNAs, and circRNAs (http://www.geneontology.org). GO analysis of differentially expressed peaks or RNAs was performed by using R software based on hypergeometric distribution. The number of genes related to altered peaks or enriched in each GO term was counted, and the related gene attributes annotated with a GO term were identified (http://www.geneontology.org).

In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (www.genome.jp/kegg) was performed by using R software, and a hypergeometric distribution test was performed to evaluate the significance of genes related to altered peaks or that of differential expression of mRNAs, lncRNAs, and circRNAs in association with each pathway term.

m6A-IP-qPCR and RT-qPCR

Twenty-two genes with differentially methylated sites as determined according to MeRIP-seq were verified by RT-qPCR. A small amount of fragmented RNA was used as the input control. The remainder of the RNA was incubated with anti-m6A antibody-coupled beads. The m6A-interacting RNAs were then immunoprecipitated and eluted from the beads. Both input control and m6A-IP samples were subjected to RT-qPCR with gene-specific primers.

Statistical analysis

All statistical analyses were performed using GraphPad Prism 8.0 software. The classification data were described as percentages, and the differences between individual groups were determined by Chi-square test. The measurement data were expressed as the mean ± standard deviation, and the significance of differences between two groups was determined by t-test or Pearson correlation test. A two-tailed P value less than 0.05 was considered to be statistically significant.