Genetic Variants in CPA6 and PRPF31 Are Associated With Variation in Response to Metformin in Individuals With Type 2 Diabetes

The incidence of type 2 diabetes (T2D) is increasing and the Centers for Disease Control and Prevention estimate that 11% of U.S. adults, aged 20 years or older, have diagnosed or undiagnosed T2D, and 35% of people in the same age-group have prediabetes based on fasting glucose or hemoglobin A_1c (HbA_1c) levels (1). T2D is now considered a global epidemic with prevalence increasing from 108 million in 1980 to 422 million in 2014 (2). Individuals with T2D are at an increased risk of developing blindness and kidney failure and are at risk for lower-limb amputations. Furthermore, individuals with T2D are 2–4 times more likely to develop cardiovascular disease (CVD), including heart attack and stroke (3).

Metformin, a member of the biguanide class of drugs, is now considered first-line therapy for individuals with T2D (4). Despite intensive investigation, the molecular mechanisms mediating metformin’s beneficial effects on glycemic control remain controversial and poorly understood. Studies performed predominantly in animal models have implicated various mechanisms and molecules as participating in metformin’s beneficial effects, including direct inhibitory effects on mitochondrial function, activation of hepatic AMPK, and alterations in glucagon signaling pathways (5–10). The relative importance of these mechanisms for metformin’s beneficial effects on glycemic control in humans is unclear.

There is significant interindividual variability in response to metformin (11–13). This suggests that a better understanding of the mechanisms by which metformin functions might allow for more tailored and precise treatment. Studies suggest that heritable factors contribute to this variability, which provides an opportunity to identify causal genetic contributors through genetic association studies (11). Recent genome-wide association studies (GWAS) have identified common variants that affected metformin response in the ATM locus, which has been shown to activate AMPK, and the SLC2A2 locus, which encodes the facilitated glucose transporter, GLUT2 (12,13). Here, we use a GWAS approach in a large cohort of individuals with T2D in the Action to Control Cardiovascular Risk in Diabetes (ACCORD) clinical trial to test for associations of both common and rare variant single nucleotide polymorphisms (SNPs) with change in HbA_1c in response to metformin treatment.

The ACCORD trial followed 10,251 participants for up to 8 years at 77 clinical centers in the U.S. and Canada to compare the benefits and risks of treatment strategies for intensively targeting glycemia, blood pressure, and dyslipidemia versus standard targets in individuals with T2D at high risk for CVD (14–17). No overall benefit and possible harms were observed in the combined primary CVD end points with intensive glucose-lowering therapy. However, significant variability in treatment response was also observed, highlighting the potential for identifying genetic markers of drug response that may lead to the development of improved and personalized treatment strategies. The results presented here point to novel mechanisms of metformin drug response and potential therapeutic targets for the treatment of T2D.

The ACCORD trial (ClinicalTrials.gov identifier: NCT00000620) was a double 2 × 2 factorial design comparing intensive versus standard treatment approaches for controlling glycemia, blood pressure, and dyslipidemia that enrolled 10,251 patients with T2D with a history of CVD or at least two known risk factors for CVD, such as documented atherosclerosis, albuminuria, dyslipidemia, hypertension, smoking, or obesity (14). Additional details about randomization and selection criteria for the various arms of the ACCORD trial can be found in the Supplementary Data. Participants in ACCORD were given an option to provide a blood sample for future genetic studies, and over 80% of participants agreed to do so. A workflow describing the selection of subjects can be found in Fig. 1. Additional demographic information can be found in Table 1.

Glycemic response to metformin was evaluated in subjects who began taking metformin while enrolled in ACCORD; subjects who reported taking metformin or another biguanide prior to enrollment in the trial were excluded. Some subjects were taking additional medications and were accounted for as described in the Supplementary Data and Graham et al. (18). Subjects were scheduled for study visits every 1 or 4 months based on randomization to the intensive or standard glycemia arms, respectively. HbA_1c was recorded every 4 months. Metformin response was calculated as on-treatment HbA_1c minus pretreatment HbA_1c. Pretreatment HbA_1c levels were recorded no more than 30 days prior to the start of metformin. On-treatment HbA_1c levels were defined as the first recorded measurement acquired after at least 90 days and no more than 270 days from the start of metformin treatment. This interval was chosen to allow for HbA_1c levels to stabilize after starting metformin, while limiting changes in glycemia levels due to additional modifications to treatment regimens. Medication compliance and other details regarding the phenotype definition can be found in the Supplementary Data.

Briefly, 6,085 unique samples from ACCORD participants who consented to genetic studies conducted by any investigator were genotyped at the University of Virginia on Illumina HumanOmniExpressExome-8 v1.0 chips (Set 1); 8,174 unique samples, including the above 6,085 samples plus 2,089 samples from ACCORD participants who consented to genetic studies only if conducted by ACCORD investigators, were genotyped at the University of North Carolina on Affymetrix Axiom Biobank1 chips (Set 2). The data were then merged, resulting in one data set consisting of 5,971 samples genotyped at a total of 1,240,656 individual SNPs, and another with an additional 2,083 samples genotyped at 583,613 SNPs. Additional details can be found in the Supplementary Data and Marvel et al. (19).

We implemented a variable selection procedure to address potential confounding. Some variables were included in the model based on previous studies or expert knowledge, whereas other variables were selected based on a backward selection approach and Bayesian information criteria to identify covariates specific to the ACCORD data set. All covariate names and descriptions before variable selection can be found in Supplementary Table 1. A substantial proportion of the cohort was taking other medications during the metformin treatment response time frame and these concomitant medications were incorporated into the model as described in the Supplementary Data. To prevent confounding due to population substructure, principal component analysis was performed based on the genotype data using EIGENSTRAT (v4.2) and is described in the Supplementary Data. The final models after covariate selection can be found in Supplementary Tables 2–4. A workflow detailing each step of the analysis can be found in Supplementary Fig. 4.

Associations between a phenotype, covariates, and single common variant (minor allele frequency [MAF] >3%) were tested using the linear regression model. Genotyped variants were tested as an additive variable using PLINK, where is the number of minor alleles for the ^th individual. Imputed variants were tested using a linear regression model in the statistical programming language, R, where is the dosage score computed from the posterior probabilities for genotypes and (20). For SNPs that were only genotyped in Set 1 subjects and were imputed in Set 2 subjects, association tests were calculated for each set separately and the results were combined by meta-analysis using PLINK (21,22). Tables and figures specify whether each result derives from a SNP that was genotyped in all subjects, imputed in all subjects, or represents a combination of genotyped and imputed data using the meta-analysis approach described.

The rare variant analysis approach has been previously described (19). Briefly, all variants (MAF ≤3%) were mapped to gene annotation from Ensemble (GRCh37.p13). A suite of five rare variant tests comprised of burden and nonburden were then used to assess associations with metformin response. We combined the set of five P values from each test into a single P value for each gene using the procedure described by Dai et al. (23). Subsequently, the combined P value was corrected for multiple comparisons with a false discovery rate approach (24,25), and a threshold of q <0.1 was used for statistical significance.

SNPs that associated with metformin response in the common variant analyses (P < 5 × 10⁻⁶) were tested for associations in the following cohorts. There was one cohort from European ancestry, which included combined samples from Genetics of Diabetes Audit and Research in Tayside Scotland (GoDARTS) and PMET1-EU (Pharmacogenomics of Metformin, cohort 1, of European ancestry; PMET1 was previously named PMT2) (REP1) (n = 6,963). There were two cohorts from African American ancestry, referred to as REP2 (n = 646) and REP3 (n = 369). Briefly, the GoDARTS cohort in REP1 has been described before by the Wellcome Trust Case Control Consortium 2 (13). The REP2 cohort was collected from two main clinical sites, Kaiser Permanente Northern California and Kaiser Permanente South East, as previously described (26). The PMET1-EU in REP1 and the REP3 cohorts were collected from the Research Program on Genes, Environment and Health (RPGEH), based at Kaiser Permanente Northern California (27). Metformin response for REP1 was defined as pretreatment HbA_1c minus posttreatment HbA_1c, resulting in effects that are in the opposite direction of those in ACCORD. Because of this, the additive inverse of effect sizes in REP1 was used in the analysis with ACCORD results. For REP2 and REP3, metformin response was calculated as posttreatment HbA_1c minus pretreatment HbA_1c.

To provide additional support for the role of the lead SNP, rs57081354, in modulating metformin action, we performed an association test between rs57081354 and baseline HbA_1c in ACCORD individuals that were already on metformin at the start of the trial (REP4) and compared those results to the association of baseline HbA_1c in individuals not on metformin at the start of the trial. All individuals in REP4 were excluded from the discovery cohort due to the lack of a recorded pretreatment HbA_1c, making it impossible to calculate a change in HbA_1c with treatment. For this analysis, an association test was performed as described in the common variant analysis with the exception of the phenotype being baseline level of HbA_1c rather than change in HbA_1c.

Each SNP with P < 5 × 10⁻⁸ in the discovery cohort or with P < 5 × 10⁻⁶ in the discovery cohort and q <0.01 in at least one replication cohort were meta-analyzed across available cohorts using PLINK (21,22). Meta-analysis with P < 0.05 was considered to be statistically significant. Additional information regarding each cohort can be found in the Supplementary Data, and a workflow describing how the replication was evaluated can be found in Supplementary Fig. 4.

Dietary questionnaires were assessed for 1,600 participants that both consented to genotyping and had dietary records as part of the Health-Related Quality of Life substudy within ACCORD. Scores in response to the dietary questionnaires were tested for association with the lead SNP, rs57081354, using a linear regression model and the same covariate selection process as described above. Details regarding this analysis can be found in the Supplementary Data.

A total of 1,312, 845, and 222 subjects were included in the common and rare variant analyses for all races combined, white, and black cohorts, respectively (Fig. 1). Other racial groups were included in the all races combined group, but sample sizes were too small to perform stratified analyses. Variation was observed in HbA_1c response (Fig. 2). The mean change in HbA_1c was −1.42% (95% CI −1.49, −1.36) for all races combined, −1.47% (−1.55, −1.39) for white subjects only, and −1.35% (−1.50, −1.19) for black subjects only. The difference in the change in HbA_1c between white and black subjects was not statistically significant (P = 0.16). Variables selected as being significantly associated with metformin response explained approximately 55% of the variation in metformin response and were subsequently included in the common and rare variant analyses (Supplementary Tables 2–4). Briefly, assignment to the intensive glycemia arm of the ACCORD trial, years with T2D, and pretreatment HbA_1c were significantly associated with metformin response in all tested groups.

A total of 316,203 genotyped, 547,639 meta-analyzed, and 7,567,403 imputed variants had MAF >3% and were included in the common variant analysis when all races were combined. A total of 1, 3, and 7 SNPs were significantly associated with HbA_1c response in all races combined, white subjects only, and black subjects only, respectively (P < 5 × 10⁻⁸). In all subjects combined, the lead SNP rs57081354 located within an intron in the NBEA gene associated with HbA_1c response (Fig. 3). Sixty-five SNPs in 16 loci on 14 different chromosomes (in genes NBEA, LOC105374140, CXCL13, MCC, FAM189A1, LOC107984525, LOC105379109, KIAA1024, RBFOX1, LOC643339, UBOX5, NR2C1, BPHL, PRPF31, LOC105377165, and ADAMTS9-AS2) reached the threshold for suggestive significance (P < 5 × 10⁻⁶). For white subjects only, 3 SNPs in KIAA1024 were associated with HbA_1c response at genome-wide levels of significance (P < 5 × 10⁻⁸), and 23 SNPs in 7 loci (in genes KIAA1024, LOC283177, LOC105377165, ADAMTS9-AS2, CPA6, NR2C1, and DPYD) reached the threshold for suggestive significance (P < 5 × 10⁻⁶) (Supplementary Fig. 1). When black subjects were analyzed separately, 2 SNPs in the LOC102724874 locus were associated with HbA_1c response at genome-wide levels (P < 5 × 10⁻⁸), and 53 SNPs in 8 loci (in genes LOC102724874, LOC105374308, LOC105369406, SLC35D3, RHPN2, ETS1, SDK2, and GPC6) reached the threshold for suggestive significance (P < 5 × 10⁻⁶) (Supplementary Fig. 2). Lead SNPs associated with HbA_1c response (P < 5 × 10⁻⁶) are presented in Tables 2–4.

Rare variants (MAF ≤0.03) in a total of 17,078 genes were tested for association with HbA_1c in all subjects combined and separately in white and black subjects. In the combined subject analysis, rare variants in STAT3 were significantly associated with HbA_1c response (q <0.1) (Supplementary Fig. 3). Five rare variants were available for testing in STAT3: rs146620441 (monomorphic), rs140604473 (missense), rs149214040 (missense), rs114401618 (intronic), and rs17882069 (synonymous) (28). Importantly, rs146620441 was monomorphic and the other four SNPs only showed variation in black subjects, although there were no statistically significant rare variant associations observed in race-based analyses (q >0.2). We tested whether STAT3 failed to yield significant associations in the stratified analysis due to the small number of black subjects (n = 222) by oversampling black subjects until we had the same number of subjects as when all races were combined (n = 1,312). Oversampling to build a larger cohort resulted in a q = 0.0016 versus a q = 1 in the original cohort of black subjects, supporting the premise that this rare variant finding is relevant for black subjects even though it was detected when all subjects were combined due to an increase in statistical power.

NBEA has been previously reported to influence sugary food preference in animal models (29), and we further assessed whether the lead SNP in this locus, rs57081354, was associated with dietary scores. Data were available in 1,553 ACCORD subjects (Supplementary Table 6). The C-allele in rs57081354 that associated with worse metformin effect also associated with increased consumption of dessert (q = 0.12). Subjects with the C-allele in rs57081354 were more likely to answer “Yes” to the question, “Did you eat dessert?” within past 3 months (β = 0.044 [95% CI 0.012, 0.075]).

A total of 134 unique SNPs were associated with metformin response (P < 5 × 10⁻⁶) when all races were combined, in white subjects only, and in black subjects only. Of these SNPs, 81, 64, and 61 were available in REP1, REP2, and REP3 cohorts, respectively. The most significant finding in the discovery cohort was rs57081354 located in the NBEA gene (P = 4.02 × 10⁻⁸, β = 0.323). This finding did not replicate in REP1 (P > 0.05) and was not available for testing in REP2 or REP3. Out of the 26 SNPs in NBEA with P < 1 × 10⁻⁶, only rs1337379 was tested in REP2 and REP3 but was not associated with metformin response (P > 0.05). SNPs in NBEA were not significantly associated with metformin response in ACCORD when black or white subjects were tested separately, only when all races were combined. Interestingly, in the set of subjects who were already on metformin when they enrolled in ACCORD (REP4), rs57081354, was marginally associated with higher baseline HbA_1c (P = 0.09), whereas rs57081354 was not associated with baseline variation in HbA_1c in those not on metformin (P = 0.63). Although the phenotypes for REP1 and REP4 are not identical, meta-analysis of the ACCORD discovery cohort with REP1 and REP4 produced an association representing higher HbA_1c on metformin (P = 0.009) (Fig. 3).

In ACCORD, rs254271, located in an intron in pre-mRNA processing factor 31 (PRPF31), was associated with worse metformin response when all races were combined (P = 3.79 × 10⁻⁶). In the REP1 cohort, rs254271 was associated with worse metformin response when corrected for baseline HbA_1c (P = 6.21 × 10⁻⁵) and without adjusting for baseline HbA_1c (P = 1.57 × 10⁻²). Genotypes for rs254271 were not available for analysis in the African American replication cohorts, REP2 and REP3. Meta-analysis of ACCORD with REP1 produced a significant association for rs254271 with worse metformin response (P = 1.2 × 10⁻⁸) (Fig. 4).

In white subjects in ACCORD, rs2162145, located in carboxypeptidase A6 (CPA6), was associated with better metformin response (P = 4.04 × 10⁻⁶). This finding was consistent with findings in African American subjects in REP2 (P = 0.006). Although, rs2162145 was not significantly associated with metformin response in black subjects in ACCORD (P = 0.14), it is possible that we were underpowered to detect an association due to the relatively small number of black subjects available (n = 222) and that meta-analysis could increase statistical power through the inclusion of additional subjects. Meta-analysis of black subjects in ACCORD and those in REP2 and REP3 demonstrated the most significant association of rs2162145 with better metformin response (P = 0.005) compared with meta-analysis of all white subjects (ACCORD white subjects and REP1) and all cohorts (ACCORD, REP1–3) with P = 0.22 and P = 0.09, respectively (Fig. 4).

T2D adversely affects the quality of life for millions of individuals and places a significant burden on the health care system in the U.S. and globally (30,31). Although many treatment options for T2D are available, metformin has remained the first-line treatment for T2D for decades. Nevertheless, the mechanisms by which metformin lowers blood glucose are not well understood, and factors that may predict which patients will have optimal responses to metformin are also not well understood. Here, we use a GWAS approach in a large cohort of individuals with T2D in the ACCORD clinical trial to test for both common and rare variant SNPs that associate with change in HbA_1c level in response to metformin treatment. To support our initial findings, we explored these associations in multiple independent cohorts.

A previous study conducted in a large cohort using electronic medical records found that African American patients on metformin had a larger reduction in HbA_1c than European Americans (32). We failed to find a significant difference between metformin response in white subjects and black subjects in ACCORD, with changes in HbA_1c of −1.47% and −1.35%, respectively (P = 0.16). It is possible that enrollment in a more structured clinical trial impacts these outcomes differently than the far less structured care reflected in electronic medical record data. Additional research will be needed to definitely characterize the response of metformin across ethnicities. Previous GWAS have identified rs11212617 in the ATM locus and rs8192675 in SLC2A2 as being significantly associated with metformin response (12,13). However, neither SNP was significantly associated with metformin response here (P > 0.1). These SNPs may have failed to replicate in the current study due to cohort differences and, in particular, the sample size of the ACCORD cohort and the use of concomitant medications. Information regarding the prescribed dose for metformin was not recorded in ACCORD and may have impacted our findings. However, as described in research design and methods, we sought a minimal time window to assess the response of metformin while limiting the likelihood of other medication changes. Additional details about how these studies compare with ACCORD can be found in the Supplementary Data.

In the analysis of subjects from all racial groups combined (n = 1,312) (Fig. 1), those with the C-allele of rs254271, an intronic variant in PRPF31, had a worse metformin response at suggestive significance levels (P = 3.79 × 10⁻⁶, β = 0.16). The effect of this variant was also observed in REP1, consisting of subjects of European ancestry (P = 6.21 × 10⁻⁵, β = −0.073). Interestingly, rs254271 is an expression quantitative trait locus in several tissues for PRPF31 and nearby genes, including NADH:ubiquinone oxidoreductase (NDUFA3), which codes for a subunit of the complex I of the respiratory chain (33–35). Metformin is known to inhibit complex I of the respiratory chain to cause AMPK activation (5). In addition, previous studies have shown that PRPF31(+/−) knockout mice have increased total body fat (P = 1.78 × 10⁻⁶) and increased fasted circulating glucose (P = 5.73 × 10⁻⁶) (36).

Genotypes for rs254271 were not available for analysis in the African American replication cohorts, REP2 and REP3. The phenotype in REP1 represents HbA_1c reduction, so the direction of the effect for rs254271 is consistent between ACCORD and REP1 and represents worse HbA_1c response in individuals with the C-allele (Fig. 4). Mutations in PRPF31 have been associated with an autosomal dominant rare Mendelian disease, retinitis pigmentosa (37). PRPF31 was previously identified as a candidate gene for diabetic retinopathy in a genome-wide linkage study in Mexican Americans with T2D (38). We tested rs254271 for association with eye disorders in ACCORD; however, no significant associations were found (P > 0.05) (data not shown). Further research will be needed to determine what role SNPs in PRPF31 may have in metformin drug response or in diabetic retinopathy.

In the analysis of only white subjects in ACCORD, rs2162145, in the CPA6 gene, was associated with metformin response at suggestive significance levels (P = 4.04 × 10⁻⁶, β = −0.197). This same SNP was associated with metformin response in the REP2 cohort (P = 0.006, β = −0.127), but not in REP1 or REP3 (P > 0.05) (Fig. 4). To our knowledge, SNPs in CPA6 have not previously been associated with metformin response; however, rs2162145 has a weak significant association with fasting glucose (P = 0.02) (39). In addition, a meta-analysis of T2D in 17 GWAS comprising 8,284 case and 15,543 control subjects in African Americans identified rs7003257, an intron variant in CPA6, as being associated with T2D (P = 1.17 × 10⁻⁶) (40), and this SNP has also been previously reported to affect 2-h glucose levels (41). We were underpowered to detect an association in black subjects. However, meta-analyzing all available black subjects in ACCORD, REP2, and REP3 improved the statistical power so an association was detectable (P = 0.005). It will be necessary for additional studies to investigate the extent to which CPA6 is involved with metformin drug response.

A genotyped SNP in the NBEA gene on chr13 rs57081354 was associated with change in HbA_1c levels at genome-wide significance (P = 4.02 × 10⁻⁸, β = 0.323) in our discovery cohort (Fig. 3). SNPs in NBEA did not replicate in REP1–3 cohorts. rs57081354 was marginally associated with baseline HbA_1c levels in individuals who were taking metformin at the time they enrolled in ACCORD, with the C-allele displaying higher HbA_1c values (P = 0.09) (Fig. 3C). NBEA is a protein kinase A anchor protein that is expressed at low levels ubiquitously and at higher levels in the brain, pituitary, and β-cells in the pancreas (42). Interestingly, haplotype-insufficient Nbea(+/−) mice prefer sugary foods (glucose and fructose) more than their wild-type controls, gain more weight per unit of food than wild-type mice, and have higher insulin and leptin levels (29). Subsequently, we tested for associations between rs57081354 and 19 dietary scores recorded for a subset of 1,553 ACCORD subjects and found that rs57081354 was associated with an increase in dessert consumption (q <0.2) in those with C-allele (Supplementary Table 6). An expanded summary of the literature regarding NBEA can be found in the Supplementary Data.

Additionally, we tested rare variants for associations with change in HbA_1c, and STAT3 when all races were combined, was the only gene that achieved statistical significance for metformin treatment response (q <0.1). However, as described in results, black subjects were the only population with rare variation in STAT3, and STAT3 was not detected in the stratified analysis likely due to the small number of subjects available for analysis (n = 222). In this same cohort, rare variants in STAT3 were not associated with baseline HbA_1c levels (q >0.2). STAT3 encodes a ubiquitously expressed transcription factor with broad activities in metabolism, immunity, and cell growth (43). STAT3 is an essential signal transducer for the key metabolic hormone leptin and the pleiotropic cytokine interleukin 6 with broad effects on feeding behavior and fuel homeostasis. Interestingly, metformin has been shown to interact with STAT3 signaling in chemotherapeutic applications and in an inflammatory bowel disease model (44–47). To our knowledge, interactions between metformin and STAT3 signaling have not been investigated with regard to glucose homeostasis.

Here, we present evidence for PRPF31, CPA6, and STAT3 being involved in novel glucose-lowering mechanisms for metformin. Additional investigation will also be required to confirm the clinical impact of the identified variants with regard to metformin therapy in subjects with diabetes. Although prior studies suggest roles for NBEA and STAT3 genes in glucose metabolism, additional research will be required to investigate mechanisms by which they might interact with metformin. Such studies are likely to provide new insight into mechanisms regulating glucose metabolism and may point the way toward novel therapeutic targets for more precise interventions in T2D.

Funding. This research was supported by National Heart, Lung, and Blood Institute grants 5R01HL110380-04 to the University of North Carolina at Chapel Hill and R01 HL110400 to the Joslin Diabetes Center and the University of Virginia. S.W.Y. and K.M.G. were supported by grants from the National Institutes of Health (GM61390 and GM117163). E.R.P. holds a Wellcome Trust Investigator Award (102820/Z/13/Z).

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. D.M.R., S.W.Y., M.A.H., H.L.M., M.J.W., J.B.B., and A.A.M.-R. wrote the article. D.M.R., M.M.H., M.K., J.C.M., H.L.M., A.D., K.M.G., E.R.P., M.J.W., J.B.B., and A.A.M.-R. designed the research. D.M.R., S.W.Y., K.Z., S.W.M., H.S.H., T.M.H., H.G., M.J.W., and A.A.M.-R. performed the research. D.M.R., S.W.Y., K.Z., S.W.M., and J.R.J. analyzed the data. A.A.M.-R. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Prior Presentation. Parts of this study were presented in abstract form at the American Society for Clinical Pharmacology & Therapeutics Annual Meeting, Washington, DC, 15–18 March 2017.