Plasma miRNA profiles associated with stable warfarin dosage in Chinese patients

Study design

We enrolled 99 patients on different stable doses of warfarin under treatment at the Cardiovascular Hospital (School of Medicine, Xiamen University). A nested case-control study was designed to identify plasma miRNAs as surrogate biomarkers for warfarin precision medicine. Bioinformatics analysis revealed the key genes related to warfarin dose, then a set of warfarin dose-related miRNAs were predicted using the miRNAs-target gene prediction databases. Our predictions were verified using stem-loop quantitative reverse transcription and hydrolysis probe-based RT-qPCR assays in 99 patients. Plasma collection and diagnosis were carried out following standard procedures. Written informed consent was obtained from all patients before study participation. The Ethics Committee of the School of Medicine, Xiamen University, China approved the protocol (IRB approval number: KY2014001), and the study was conducted following the principles of the Declaration of Helsinki.

Sample collection and patient characteristics

Patients were eligible if they were 18 years of age or older and on long-term warfarin therapy. Other eligibility criteria included: (1) Chinese Han ethnicity, (2) baseline liver and kidney function test data and international normalized ratio (INR) results available, (3) not taking other anticoagulant drugs, (4) not undergoing chemotherapy, (5) no family history of alcoholism, and (6) no history of chronic liver failure or other liver diseases.

According to the consensus of domestic experts (Chinese Society of Cardiology, 2013), the standard targeted INR is 2.0–3.0, but it is 2.5–3.5 for patients with acute valve replacement and two-valve implantation. A stable dose is defined as the mean dose of three consecutive measurements taken at intervals of more than one week to determine the same daily dose to maintain INR within the target range. The study subjects were divided into low-dose (≤2 mg/day), medium-dose (2 to 4 mg/day), and high-dose (≥4 mg/day) groups. Among the 99 patients, 30 (30.3%), 40 (40.4%), 29 (29.3%) were classified as low-, medium-, and high-dose, respectively. patient demographic and clinical characteristics were obtained from the Laboratory Information System or a follow-up telephone call (Table 1).

Bioinformatics analysis

PharmGKB is a pharmacogenetics and pharmacogenomics database detailing how genetic variation affects drug response. It provides a theoretical basis for research into personalized medicine and pharmacogenomics based on the assumption that differential expression of key pathway genes may impact dose requirements. We screened out key genes that were supported at least two annotations in PharmGKB.

For functional enrichment analysis of predicted genes, GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analyses were conducted using the DAVID (Database for Annotation, Visualization, and Integration Discovery, https://david.ncifcrf.gov/summary.jsp). Parameters for the enrichment analysis were as follows: p < 0.05, gene count ≥2, top 5.

Five miRNAs prediction sites (miRanda, miRDB, miRWalk, RNA22, TargetScan) were used to identify relevant miRNAs that target warfarin key genes. The final result set was validated in subsequent RT-qPCR assays. The commonly used miRNAs prediction sites were summarized in Table 2.

miRNAs isolation and quantification

All patients in this study had peripheral blood drawn at study enrollment. Blood samples were collected in commercial vacuum blood collection tubes containing 3.2% sodium citrate and centrifuged at 3500 rpm for 10 min to obtain plasma. After 15-min high-speed centrifugation at 13300 rpm at room temperature to remove cell debris, the supernatant plasma was recovered and stored at −80 °C until analysis.

To assess the expression level of the targets and minimize operational bias, cel-miR-39 was used as an external reference. Before miRNAs extraction, 13.2 pg cel-miR-39 was added to each plasma sample, followed by the same extraction and amplification process as the targets. MiReasy Serum/Plasma kits (cat. #217184; Qiagen, Hilden, Germany) were used to extract miRNAs from 200 µl plasma (including 13.2 pg cel-miR-39) with column extraction, then 25 µl elution buffer was added. Next, we carried out TaqMan probe-based RT-qPCR assays with two commercial kits on the Applied Biosystems 7500 Sequence Detection System (Applied Biosystems, Foster City, CA, USA): TaqMan™ MicroRNA Reverse Transcription Kit based on the stem–loop primer method (cat. #4366597; Applied Biosystems) for the reverse-transcription step and TaqMan universal master MixII, no UNG (cat. #4440040; Applied Biosystems) for qPCR assays. All reactions, including the no-template controls, were run in duplicate. All operations were carried out following standard procedures.

Statistical analyses

All statistical analyses and plots were conducted using IBM SPSS Statistics 25 (IBM Corp, Armonk, NY, USA). Baseline characteristics were compared by the Kruskal–Wallis test (for continuous variables) or Fisher’s exact test (for categorical variables) to analyze differences among the three groups. The relative expression level of target miRNAs (denoted by ΔΔCt) was corrected by cel-miR-39 (ΔCt _sample= Ct_target −Ct_{cel−miR−39}) and the reference sample (ΔΔCt _sample = ΔCt_sample − ΔCt_{ref. sample}) and then used as the data for subsequent statistical analyses. We used Kruskal–Wallis test and multivariate logistic regression to compare differences in plasma miRNAs relative expression level among the three warfarin dosage groups. For all tests, two-tailed p < 0.05 was considered statistically significant.

Patient characteristics

From April 2018 through November 2018, we enrolled a total of 99 subjects. We recorded the baseline demographics and clinical characteristics of the patients receiving stable warfarin treatment. There were no significant differences in baseline characteristics among the low-, medium- and high-dose groups (Table 1).

Identification of key genes

Overall, 52 relevant genes were identified in the PharmGKB database, of which 31 genes were supported by at least two annotations (Table 3). Three or more evidences support the CYP2C9, VKORC1, CALU, CYP2C19, CYP1A1, CYP4F2, EPHX1, GGCX, and PROS1 genes. PharmGKB clinical annotations (CAs) provide information about variant-drug pairs based on variant annotations in the database. Scientific curators manually review variant annotations and create genotype-based summaries describing the variants’ phenotypic impacts (Whirl-Carrillo et al., 2012). CAs were divided into four levels according to the strength of the evidence, with CA¹ meeting the highest criteria. VKORC1, CYP2C9, and CYP4F2 genes were rated as CA¹, and their effects on warfarin dosage have been widely reported (Tavares, Marcatto & Santos, 2018).

To better understand the functions of the 31 genes predicted by PharmGKB, we performed GO enrichment and KEGG pathway analyses. In the GO enrichment analysis, the predicted genes were significantly enriched in biological processes related to blood coagulation, hemostasis, and body fluid regulation. KEGG pathway analysis revealed that the predicted genes were mainly enriched in chemical carcinogenesis, retinol metabolism, and complement and coagulation cascades (Fig. 1).

According to the PharmGKB “PW” annotation and the results of pathway enrichment analyses, we divided the 31 genes into two groups: pharmacokinetics and pharmacodynamics. We followed conditions to select the key genes of each pathway for subsequent miRNAs prediction. The conditions were as follows: (1) at least two supporting annotations, and (2) top 2 levels for “CA” or the widely reported genes related to warfarin dosage. Finally, we obtained three key genes related to pharmacokinetics: CYP1A1 (Luo et al., 2017), CYP2C9, and CYP2C19 (Zhang et al., 2016), and four key genes in pharmacodynamics: VKORC1, GGCX, CYP4F2, and CALU.

Prediction of miRNAs targeting key genes

Commonly used miRNAs prediction databases include miRanda, miRDB, miRWalk, RNA22, TargetScan. MiRWalk provides more comprehensive predictions because it integrates two predicted databases (TargetScan and miRDB) and one validated database (miRTarBase). We used miRWalk as the primary tool, with the others as supplements. We inputted the 31 differential genes and screened out miRNAs with p < 0.05. We followed three principles for further screening: (1) binding with 3′ UTR, CDS or 5′ UTR, the more binding sites, the higher score the miRNA was; (2) consistent output in at least two databases and. (3) miRNAs that can target at least three key genes of warfarin. They were summarized in Table 4. The key genes related to drug metabolism and drug effects obtained from the above results were used for further prediction (Fig. 2).

About warfarin pharmacodynamics, we obtained five of the most relevant miRNAs: hsa-miR-937-5p, hsa-miR-4700-5p, hsa-miR-1276, hsa-miR-500a-3p, and hsa-miR-940. For warfarin pharmacokinetics, hsa-miR-211-3p and hsa-miR-6515-5p were the most relevant. Warfarin exerts its anticoagulant effect as an inhibitor of VKOR, and it is metabolized in the liver by the enzymes encoded by the cytochrome P450 (CYP) superfamily of genes. The CYP superfamily is responsible for >75% of phase drug metabolism in the liver, so other prescriptions taken by enrolled patients may influence warfarin pharmacokinetics-related miRNAs. Moreover, as a direct target of warfarin, VKORC1 plays a leading role in individualized warfarin dosage. We therefore selected VKORC1, which highly related to warfarin dose, for further screening and verification. We inputted “VKORC1” and screened out miRNAs with p < 0.05. miRNAs predicted by at least two databases were considered as VKORC1-related miRNAs. Next, we used miRWalk to reversely predict the targeted genes of these miRNAs, finally screening out 20 miRNAs related to VKORC1 (Table 5). Notably, hsa-miR-1276 was the only one with more than one binding site on VKORC1 mRNA. Synthesizing the above results and information in the literature (e.g., essential functions reported, expression level, etc.), we finally screened out the top three highest score miRNAs: hsa-miR-24-3p(MIMAT0000080), hsa-miR-133b(MIMAT0000770), and hsa-miR-1276(MIMAT0005930).

Validation study

We examined the three predicted miRNAs by hydrolysis probe-based RT-qPCR in the three dosage groups. Kruskal–Wallis test showed that hsa-miR-133b (p < 0.001) and hsa-miR-1276 (p = 0.024) were significantly different among the three groups, but hsa-miR-24-3p (p = 0.806) was not.

Logistic regression analysis was used to adjust for confounding factors. After adjusting for clinical covariates, including age, gender, height, weight, smoking habit, alcohol intake, and combination use of amiodarone, there were significant between-group differences among the three-dose groups for hsa-miR-133b expression (p = 0.005). However we observed an “n-shaped” dose-dependent curve rather than a linear relationship. There were no significant differences in hsa-miR-24-3p (p = 0.475) or hsa-miR-1276 (p = 0.558) among the three groups (Fig. 3).