Methylation of CYP1A1 and VKORC1 promoter associated with stable dosage of warfarin in Chinese patients

Study subjects and study design

We undertook a two-stage case–control study that enrolled patients who were on stable warfarin treatment from a biobank(Planning and preparation in 2010, collection between 2011 and 2018) of Chinese patients treated at the Zhongshan Hospital, School of Medicine, Xiamen University and the Cardiovascular Hospital, School of Medicine, Xiamen University. The Ethics Committee of the School of Medicine, Xiamen University, China approved this study protocol(IRB approval number: KY2014001), and the study was conducted in accordance with the principles of the Declaration of Helsinki.

Patients on stable warfarin treatment were eligible for inclusion in the study if they fulfilled the following criteria: (1) of Chinese ethnicity; (2) requiring long-term warfarin anticoagulant therapy (target INR 2.0–3.5); (3) were 18 years of age or older; and (4) provided written informed consent. Exclusion criteria were: (1) chronic liver failure; (2) use of other anticoagulant drugs; (3) other drugs known to interact with warfarin (according to the report of Gage et al. (2008), eg. Losartan) with the exception of amiodarone, (4) currently undergoing chemotherapy; and (5) alcoholism.

We enrolled 161 patients during the study period. At the first stage, we included 30 patients as the “Screening cohort” to measure genome-wide DNA methylation and screen for potential warfarin dose-associated CpGs among the three groups (low-dose: 10, medium-dose: 10; high-dose: 10). Additionally, at the second stage, 131 patients were enrolled as the “Validation cohort” in the subsequent pyrosequencing to further study the association between the stable warfarin dose and CpG methylation among the three groups (low-dose: 37; medium-dose: 59; high-dose: 35). The study design and patient disposition is shown in File S1.

Data collection and patient characteristics

Baseline demographics (age, sex, height, weight, and habits of smoking and alcohol intake) and clinical characteristics (combination use of amiodarone, indications for warfarin treatment, and stable dose of warfarin) of study subjects (Table 1) were obtained from electronic medical records and through a follow-up telephone call.

Patients with INR in the target range had a minimum interval of 4 weeks between INR measurements. The maximum interval between INR measurements was 12 weeks, and was seen in patients on VKA therapy with consistently stable INRs. Stable warfarin dosage was defined on the basis of three consecutive INR measures in the therapeutic range (2.0–3.5) without any adjustment of warfarin dosing between the minimum and maximum intervals. The mean dose recorded during consecutive visits was defined as the stable dose of warfarin.These definitions were chosen refer to the previous studies (Cai et al., 2017; Tavares et al., 2018).

The study subjects were assigned to “low-dose” (≤2 mg/day), “medium-dose” (2 mg/day to 4 mg/day), and “high-dose” (≥4 mg/day) groups, on the basis of the stable dose of warfarin refer to the previous study in China (Cai et al., 2017). The low dose group and high dose group were defined as case groups,and the medium dose group was defined as control group.

Sample collection

All patients in this study underwent peripheral blood sampling at the time of study enrollment, which were derived from the remaining samples of medical examination. Blood specimens in commercial vacuum blood collection tubes containing 3.2% sodium citrate were separated by centrifugation at 3500 rpm for 10 min to acquire the peripheral blood cells. After centrifugation, the acquired specimens were stored at 4 degrees freezer. Then, the acquired specimens were transferred to -80 degrees ultra-low temperature freezer in 3 h. The specimens were stored in the ice box during the transfer process. All of the steps(from sampling to freeze) were completed in 6 h.

Genome-wide DNA methylation analysis

From the screening cohort, on the basis of the stable warfarin dose, three groups (low-dose: 10, medium-dose: 10;high-dose: 10) were matched by the patients, characteristics. We screened for potential warfarin dose-associated CpGs by comparison of genome-wide DNA methylation between the case groups(low dose group and high dose group) and control groups(medium dose group).

We isolated DNA from peripheral blood cell samples by using the QIAamp DNA Blood Mini Kit (Cat.#51306; Qiagen, German). DNA bisulfite conversion was done in a bisulfite batch with the Zymo EZ DNA MethylationTM Kit (Cat.#D5001, Zymo, USA). Genome-wide DNA methylation was measured by Illumina Infinium HumanMethylation 450 K BeadChip (Cat.#WG-314-1001, Illumina, USA).

Data from genome-wide methylation were analyzed using the GenomeStudio Software (version V2011.1, Illumina, San Diego, CA, USA). The average β-value was estimated as the proportion of total signal intensity from the methylated-specific probe to represent the methylation level (range 0–1). Delta beta (Δβ) was used to represent the difference in methylation levels among the case groups(low dose group and high dose group) and control group(medium dose group). To avoid false positives, the probes were filtered out using a detection P-value of less than 0.05. Moreover, probes on the X and Y chromosome were removed to eliminate sex bias. The FDR-adjusted DiffScores were computed to account for multiple testing as well as to limit false positives. The differences in methylation levels among the three dose groups were tested using the Illumina Custom model in GenomeStudio software (Carless, 2015). Significant differences were established with a FDR adjusted DiffScore $\ge \left|13\right|\tilde{P}\le 0.05$. The genome-wide methylation analysis was conducted refer to the previous pipelines and workflows (Maksimovic, Phipson & Oshlack, 2017; Morris & Beck, 2015; Morris et al., 2014).

The heatmap of hierarchical clustering for samples and CpGs (FDR adjusted DiffScore $\ge \left|33\right|\tilde{P}\le 0.001$) was generated by the Ward method using the R package pheatmap (version 1.0.12). The Volcano plot depicts the genome-wide distribution of hypo- and hypermethylated CpGs based on their Δ β and DiffScore using the R package ggplot2 (version 2.2.1). The bar chart represents the distribution of the differential CpGs.

Gene ontology, pathway, and drug database enrichment analysis

The genes, including differential CpGs, were uploaded to the Web-based gene set analysis toolkit (WebGestalt) to carry out the gene ontology (GO), the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database, and the drug database (DrugBank and GLAD4U) enrichment analysis (Liao et al., 2019). The enrichment analysis was conducted based on the parameters(minimum number of IDs in the category: 5; maximum number of IDs in the category: 2000; FDR Method: BH; significance level: Top 10).

To further undertake dimension reduction and screen for the differential CpGs, we used a rapid review (PubMed: “warfarin” AND (“gene” OR “SNP”)) to collect genes that have been shown to encode proteins of the warfarin-associated pharmacokinetics and pharmacodynamics pathways. Then, we obtained the intersection of differential CpGs and warfarin-associated genes, and CpGs within the intersection were chosen as target CpGs for further validation. The gene list is shown in File S1.

Pyrosequencing assay

To verify some of the significant CpGs from the genome-wide analysis, we selected the target CpGs at the warfarin-associated genes. Thereafter, we used pyrosequencing to detect the methylation levels of CpGs with more samples in the validation cohort (low-dose group: 37; medium-dose group: 59; and high-dose group: 35).

DNA was extracted from peripheral blood cells using the TIANamp Blood DNA Kit (Cat. #DP318, Tiangen, China) and DNA bisulfite conversion was conducted in a bisulfite batch with the EpiTect Bisulfite Kit (Cat. #59104, Qiagen, German); then, the converted DNA used as the template for PCR amplification. All primers for PCR and sequencing were designed by PyroMark Assay Design 2.0 Software (version 2.0.01.15; Qiagen, German). Pyrosequencing was carried out by PyroMark Q24 (Qiagen, German) (Clark et al., 2006; Tost & Gut, 2007). The EpiTect Control DNA and Control DNA Set(Cat. #59695, QIAGEN; German) was used in the pyrosequencing. In addition,the PyroMark Q24 Validation Oligo(Cat. #979203, QIAGEN; German) was also used to check the performance of the PyroMark Q24 system. The protocol of PCR and primers is described in File S1.

Statistical analysis

Normality of distribution was assessed with the Shapiro–Wilk normality test. The uncertain distribution of data was also considered. Continuous variables were expressed as mean ± standard deviation (SD) or median and interquartile range (IQR) when appropriate. Categorical variables were expressed as numbers and percentages.

The sample size of low dose group and high dose group was less than 50, and the distribution of the continuous variables (eg. CpGs methylation) was uncertain. Moreover, different log transformations (log10, log2, and ln) were used in warfarin stable dose to observe whether this variable transforms to the normality, but this variable was still non-normally distributed. In this case, parametric tests may lead to wrong conclusion. Therefore, the non-parametric tests were used in this study after careful consideration.

Kruskal–Wallis test (for continuous variables) and Fisher’s Exact test (for categorical variables) were used to analyze the differences between the groups (high-, medium-, and low-dose). To adjust for the effect of clinical covariates, including age, sex, height, weight, habits of smoking and alcohol intake, and combination use of amiodarone, to prevent any spurious associations, we used ordinal logistic regression analysis to evaluate the association of every CpG methylation level with the warfarin dose groups using the R package rms (version 6.1.0). Spearman rank correlation analysis was used for the correlation between every CpG methylation level and stable warfarin dose, in differential dose groups, respectively. Using the R package ggplot2 (version 2.2.1), the violin plot and dot plot were combined as a violin–dot plot to represent the differences and distributions of the CpG methylation levels among the three dose groups. The correlation map visualized the results of correlation analysis by using the R package corrplot (version 0.84).

A two-tailed P-value less than 0.05 was considered statistically significant in all tests. All statistical analyses (except the genome-wide methylation difference analysis) and plots were conducted using R (Version 3.4.4, R Foundation for Statistical Computing, Vienna, Austria).

Patient characteristics

From January 2010 through November 2018, we assigned a total of 161 patients in this study. The baseline demographics and clinical characteristics of the patients included in the screening (n = 30) and validation (n = 131) cohorts are described in Tables 1 and 2. There were no significant between-group differences in baseline characteristics in the two cohorts.

Genome-wide DNA methylation analysis

On the array, 473,920 of the 485,577 probes (97.60%) were detected with a detection P-value of 0.05 or less, after the probes on the X and Y chromosome were removed. Using differential methylation analysis, we found that 42,549, 26,076, and 3955 CpGs were differentially methylated among the three comparisons (low-dose group−high-dose group; medium-dose group−high-dose group; low-dose group–medium-dose group; FDR-adjusted DiffScore $\ge \left|13\right|\tilde{P}\le 0.05$). From the perspective of the distribution of differential CpGs on different gene regions, the differential CpGs among the three comparison pairs are mainly distributed on the Body and IGR regions. In addition, the 1stExon and the 3′ UTR regions are the regions with the least distribution of differential CpGs. From the perspective of differential CpGs on CpG island areas, the differential CpGs among the three comparison pairs are mainly distributed on the opensea area. From the perspective of the distribution of differential CpGs on chromosomes, the differential CpGs among the three comparison pairs are mainly distributed on chromosomes 1 and 6, but less on chromosomes 18 and 21. The detailed distributions of hyper- and hypomethylated CpGs among all comparisons are shown in File S2.

Furthermore, the selected differential CpGs (FDR-adjusted DiffScore $\ge \left|33\right|\tilde{P}\le 0.001$) could better differentiate among each comparison through a hierarchical cluster approach (Fig. 1A). The genome-wide distribution of hypo- and hypermethylated CpGs is shown in Fig. 1B. As with the distribution of the genomic region, when we compared the low-dose group with the high-dose group, 7950 differential CpGs, which contained 4098 hyper- and 3852 hypomethylated CpGs, were enriched in the promoter region (Fig. 1C). Most of differential hypermethylated CpGs (54.38%; 16,588/30,506) were located in the open sea region, whereas most of the differential hypomethylated CpGs (40.96%; 4934/12,043) were located in the island region (Fig. 1D).

Gene ontology, pathway, and drug database enrichment analysis

To better understand the function of the methylation-altered genes revealed in this study, we conducted enrichment analyses of the GO and KEGG pathway database and the drug database (DrugBank and GLAD4U).

In the GO enrichment analysis, we found 10 biological process terms were significantly enriched (FDR-adjusted P-value <0.05) among the 8678 differential genes associated with a warfarin stable dose. All the enrichment analyses are shown in File S3.

To further undertake dimension reduction and screening of the differential CpGs, we undertook a rapid review and verification of the results of genome-wide DNA methylation and chose cg13570656, cg12101586, and cg22549041 –three of the significantly differential CpGs at the CYP1A1 (Cytochrome P450, Family 1, Member A1) gene and cg04263740 at the VKORC1L1 (vitamin K epoxide reductase complex, subunit 1-like 1) gene as the target CpGs. CYP1A1 has been investigated to ascertain its role in R-warfarin metabolism (Zhang et al., 1995). Moreover, human VKORC1L1 was confirmed to be involved in the vitamin K metabolism system, where it can reduce inactive vitamin K 2,3-epoxide to active vitamin K; this system is the target of warfarin that is widely used in thrombolytic therapy and prophylaxis (Czogalla et al., 2018).

In an exploratory attempt in this validation study, as a post hoc analysis, we considered cg09155044 at the promoter of the VKORC1 (vitamin K epoxide reductase complex subunit 1) gene. There are three reasons why we chosen cg09155044 and the surrounding CpGs in this study, although there is no significant difference between the groups on genome-wide DNA methylation analysis. Firstly, DNA methylation at the promoter has been shown to directly silences gene in previous study (Jones, 2012). Secondly, previous study (Oldenburg, Watzka & Bevans, 2015) have shown that warfarin reduces the available reduced vitamin K and reduces the synthesis of active coagulation factors by inhibiting its target VKORC1. Finally, the cg09155044 and the surrounding CpGs are all located on the promoter of the VKORC1. Based on these considerations, we assume that methylation levels of cg09155044 and the surrounding CpGs may be related to warfarin. To verify our hypothesis, we chosen these CpGs in the validation study as a post hoc analysis. Table 3 presents detailed information on the target CpGs.

Validation study

To verify some of the significant differentially methylated CpGs from genome-wide analysis, we used pyrosequencing technology to test methylation levels of the target CpGs, from more samples in the validation cohort. In addition, the surrounding CpGs (CA-CB, CF-CG, VA, and VC-VE) were taken into the assay range of pyrosequencing (Table 4 presents detailed information on CpGs).

Table 4 shows results of the Kruskal–Wallis test and ordinal logistic regression analysis. The DNA methylation levels of CpGs (CA-CC, CE, and CG) were significantly different between the differential dose groups (P < 0.05); methylation levels of these CpGs at the CYP1A1 gene increased with an increase in the stable warfarin dose.

After adjusting for clinical covariates, there were significant between-group differences among the three dose groups for several CpGs (CA, CB, VD and VE) at the CYP1A1 gene and VKORC1 gene. However, the CpG(VKL) at the VKORC1L1 gene showed no significant difference among the three dose groups in the Kruskal–Wallis test and on logistic regression analysis.

Differences of CpG methylation levels (CA-CG and VA-VE) among the three dose groups are shown in Fig. 2A. Table 5 shows the results of Spearman rank correlation analysis. The methylation levels in several CpGs (CA-CB, CD-CE, and CG) at the CYP1A1 gene presented a positive correlation (r > 0, P < 0.05) among samples from subjects on the stable warfarin dose. Furthermore, the CpGs (CA-CE, and CG) presented a positive correlation (r > 0, P < 0.05) between the methylation levels and different dose groups. At the VKORC1 gene, the methylation levels in VD presented a negative correlation (r < 0, P < 0.05) among the different dose groups. However, the CpG (VKL) methylation at VKORC1L1 did not correlate (P > 0.05) with the stable warfarin dose variation.The correlation of the CpG (CA-CG and VA-VE) methylation levels among the three dose groups are shown in Fig. 2B.