Common genetic variation in the glucokinase gene (GCK) is associated with type 2 diabetes and rates of carbohydrate oxidation and energy expenditure

Aims/hypothesis Glucokinase (GCK) plays a role in glucose metabolism and glucose-stimulated insulin secretion. Rare mutations in GCK cause MODY. We investigated whether common variation (minor allele frequency ≥0.01) in GCK is associated with metabolic traits and type 2 diabetes. Methods Four exonic single-nucleotide polymorphisms (SNPs) and three SNPs predicted to cause loss of promoter function were identified in whole-genome sequence data from 234 Pima Indians. These seven tag SNPs and rs4607517, a type 2 diabetes variant established in other studies, were analysed in 415 full-heritage non-diabetic Pima Indians characterised for metabolic traits, and 7,667 American Indians who had data on type 2 diabetes and BMI. Results A novel 3′ untranslated region (3′UTR) SNP, chr7:44184184-G/A, was associated with the rate of carbohydrate oxidation post-absorptively (β = 0.22 mg [kg estimated metabolic body size (EMBS)]−1 min−1, p = 0.005) and during a hyperinsulinaemic–euglycaemic clamp (β = 0.24 mg [kg EMBS]−1 min−1, p = 0.0002), the rate of carbohydrate oxidation in a respiratory chamber (β = 311 kJ/day, p = 0.03) and 24 h energy expenditure, which was attributable to the thermic effect of food (β = 520 kJ/day, p = 3.39 × 10−6). This 3′UTR SNP was also associated with diabetes (OR 1.36, 95% CI 1.11, 1.65, p = 0.002), where the A allele (allele frequency 0.05) was associated with a lower rate of carbohydrate oxidation, lower 24 h energy expenditure and higher risk for diabetes. In a Cox proportional hazards model, a rate of insulin-stimulated carbohydrate oxidation lower than the mean rate at baseline predicted a higher risk for developing diabetes than for those above the mean (hazard rate ratio 2.2, 95% CI 1.3, 3.6, p = 0.002). Conclusions/interpretation Common variation in GCK influences the rate of carbohydrate oxidation, 24 h energy expenditure and diabetes risk in Pima Indians. Electronic supplementary material The online version of this article (doi:10.1007/s00125-014-3234-8) contains peer-reviewed but unedited supplementary material, which is available to authorised users.


Introduction
Glucokinase (GCK) is a hexokinase isozyme (hexokinase IV) that catalyses glucose to glucose-6-phosphate (G6P) and is involved in the first step of both glycolysis and glycogen synthesis. GCK is predominantly expressed in hepatocytes and pancreatic beta cells, with isoforms distinct in the N terminus. The pancreatic beta cell isoform is a key enzyme in regulating glucose-stimulated insulin secretion and is considered to be a glucose sensor. The liver isoform plays a central role in regulating glucose homeostasis and is a major component of the hepatic glucose-sensing system involved in glucose synthesis, breakdown and storage [1][2][3]. Rare heterozygous inactivating mutations in GCK cause MODY, mainly due to a reduced glucose-stimulated insulin secretion [4]. While rare mutations in GCK cause MODY, common variants have been associated with HbA 1c levels, fasting glucose concentrations and type 2 diabetes in white and other populations [5][6][7]. No rare coding variants in GCK were identified in 234 Pima Indians with whole-genome sequence data (unpublished data, Y. L. Muller). Thus, in the current study, we investigated the effects of common and low-frequency GCK variants with a minor allele frequency (mAF) ≥0.01 on metabolic traits and type 2 diabetes risk in Pima Indians.

Methods
Participants with outpatient longitudinal data on type 2 diabetes and BMI Electronic supplementary material (ESM) Fig. 1 shows a flow chart depicting the study design and selection of participants. All individuals in this study are participants of a longitudinal study of the aetiology of type 2 diabetes among the Gila River Indian Community in Arizona, where most of the residents are Pima Indians or Tohono O'odham (a closely related tribe) [8]. Diabetes was determined by prior clinical diagnosis or an oral glucose tolerance test according to the criteria of the American Diabetes Association [9]. A population-based sample of full-heritage Pima Indians (n=3,604, including 736 sibships [sibship is defined as sibs ≥2], Table 1) was initially used to assess associations with type 2 diabetes. A non-overlapping sample of mixed-heritage American Indians from the same longitudinal study (n=4,063, including 739 sibships; reported heritage, on average, was one-half Pima and three-quarters American Indian, Table 1) was used to assess replication. Among these samples, BMI was measured at biennial examinations and maximum BMI observed in the longitudinal study was analysed in 3,391 full-heritage Pima Indians and 3,406 mixed-heritage American Indians (Table 1) who were examined when aged ≥15 years. Fasting serum glucose concentrations were measured in 2,542 full-heritage Pima Indians and 2,887 mixed-heritage American Indians that were non-diabetic, including individuals who subsequently developed diabetes and those who remained non-diabetic ( Table 1).
Subset of participants with additional inpatient data on quantifiable metabolic traits Among the full-heritage Pima Indians described above, 415 non-diabetic individuals (including 99 sibships; male sex 58%, age 27±6 years and BMI 34±8 kg/m 2 at the time of metabolic testing) had undergone detailed studies of metabolic and anthropometric phenotypes for risk factors related to type 2 diabetes and obesity. Body composition, including percentage body fat, fat mass and fat-free mass, was estimated by underwater weighing until 1996 and by dual energy x-ray absorptiometry (DPX-1; Lunar Radiation Corp., Madison, WI, USA ) thereafter [10]. Glucose tolerance was determined by a 75 g OGTT, with measurements of fasting, 30, 60, 120 and 180 min plasma glucose and insulin concentrations [11]. A hyperinsulinaemic-euglycaemic clamp (insulin infusion rate of 40 mU m −2 min −1 with simultaneous glucose tracers) was used to measure rates of post-absorptive (basal) and insulin-stimulated glucose disappearance as previously described [11]. Indirect calorimetry measurements using a ventilated hood system were performed before and during the insulin infusion to assess rates of energy expenditure and substrate oxidation [12,13]. Pancreatic beta cell function was assessed by the acute insulin response (AIR) after a 25 g intravenous glucose bolus and calculated as the mean increment in plasma insulin concentrations from 3 to 5 min [11].
To measure 24 h energy expenditure, study participants entered a respiratory chamber for 23 h and 15 min after an overnight fast and after at least 3 days of a weight-maintaining diet [14]. Four meals were provided at 08:00, 11:00, 16:00 and 19:00 hours. Fresh air was drawn through the chamber, and CO 2 production and O 2 consumption were measured and calculated every 15 min and extrapolated to the 24 h period [15]. Spontaneous physical activity (SPA) was detected by radar sensors and expressed as percentage of time in motion per 15 min interval. The energy cost of SPA was calculated as the product of average SPA over 24 h and the slope of the regression line between energy expenditure and SPA between 08:00 and 23:00 hours [15]. Sleeping metabolic rate was defined as the average energy expenditure of all 15 min periods between 1:00 and 5:00 hours during which SPA was <1.5%, and was extrapolated to 24 h [16]. The  Allelic specific gene expression Total RNA was isolated from percutaneous abdominal adipose tissue biopsies from 14 individuals who were heterozygous for the 3′ untranslated region (3′UTR) SNP chr7:44184184-G/A. RNA was reverse-transcribed to cDNA and allelic specific expression was performed on an ABI Prism 7900 using a TaqMan probe (Applied Biosystems). The ratio of allele G/A expression was normalised to the genomic DNA control.
Statistical analyses Statistical analyses were performed using the software of the SAS Institute (version 9.2, Cary, NC, USA). A logistic regression analysis was used to assess the association of genotypes with type 2 diabetes and was adjusted for age, sex, birth year and heritage as covariates. The model was fit with the generalised estimating equation (GEE) to account for dependence among siblings. Genotype was analysed as a numeric variable representing 0, 1 or 2 copies of a given allele. To estimate the proportion of European ancestry, 45 informative markers with large differences in allele frequency between populations [17] were used as a covariate in these analyses. The association of quantifiable traits with genotypes was analysed by linear regression using the GEE procedure to account for correlation among siblings. Results were adjusted for covariates as indicated. The rate of carbohydrate oxidation, as a metabolic predictor of diabetes, was assessed in individuals with normal glucose tolerance (NGT) who had measures of carbohydrate oxidation rate at baseline during a hyperinsulinaemic-euglycaemic clamp. These individuals also had biennial follow-up OGTTs to determine diabetes status. The Cox proportional hazards model was used to determine the hazard rate ratio (HRR) for developing type 2 diabetes associated with the rate of insulin-stimulated carbohydrate oxidation including age, sex, percentage body fat, AIR and non-oxidative glucose disposal rate as covariates. Follow-up time was defined as the time from the measure of carbohydrate oxidation rate during an insulin clamp to either type 2 diabetes onset or the last evaluation when an individual remained non-diabetic. To analyse the effect of genotype on the energy expenditure trajectory, a mixed model analysis was used including time, time 2 and time 3 as fixed effects to model the non-linearity of the trajectory. Results were adjusted for age, sex, fat mass, fat-free mass and SPA.

Results
Association of GCK SNPs with carbohydrate oxidation rate Whole-genome sequence data from 234 Pima Indians were used to identify common SNPs (mAF≥0.01). Five SNPs were identified in exon regions: two synonymous amino acid substitutions, GCK-G193G and GCK-Y215Y (rs144723656), and three 3′-UTR SNPs, rs13306388, rs55714218 and a novel SNP chr7:44184184-G/A. Four SNPs were identified in the putative promoter region and were predicted to cause loss of function by the Ingenuity Variant Analysis (https://variants. ingenuity.com): rs1799831, rs1799884, rs193226243 and rs1476891. These nine SNPs and rs4607517 (which has been associated with type 2 diabetes in other studies [5][6][7] and maps 6.6 kb upstream of GCK) were analysed for pairwise LD. Eight tag SNPs capture all ten SNPs (r 2 ≥0.8, ESM Fig. 2). These eight tag SNPs (mAF≥0.01) were genotyped in a population-based sample of 3,604 full-heritage Pima Indians of which 415 non-diabetic individuals had detailed measures of metabolic and anthropometric phenotypes for risk factors related to type 2 diabetes and obesity. Associations of these eight tag SNPs with metabolic traits were analysed in these 415 individuals. For the novel 3′UTR SNP chr7:44184184-G/ A (mAF 0.05), only two of the 415 individuals were homozygous for the minor A allele, so their data were combined with those of the G/A heterozygotes for statistical analyses. Individuals with the A allele had a lower mean rate of basal carbohydrate oxidation (Fig. 1, β=0.22 mg [kg estimated metabolic body size (EMBS)] −1 min −1 per risk allele, p= 0.005, adjusted for age, sex and percentage body fat) and lower rate of insulin-stimulated carbohydrate oxidation (β= 0.24 mg [kg EMBS] −1 min −1 , adjusted p=0.0002) when compared with individuals with the G allele. However, the nonoxidative glucose disposal rate at baseline and during insulin stimulation was not different between the two groups (adjusted p=0.18 at basal state, p=0.92 during insulin stimulation).
Furthermore, individuals with the A allele had a higher lipid oxidation rate at baseline by 0.08 mg [kg EMBS] −1 min −1 ( Table 2, p=0.007, adjusted for age, sex and percentage body fat) and during insulin infusion by 0.09 mg [kg EMBS] −1 min −1 (adjusted p=0.01) when compared with individuals with the G allele. These changes in lipid oxidation rate may be secondary to the changes in carbohydrate oxidation rate. Individuals with the A allele also had a lower basal rate of endogenous glucose production (β=0.08 mg [kg EMBS] −1 min −1 , adjusted p=0.03), but no difference in endogenous glucose production rate during insulin infusion (adjusted p=0.31). This 3′UTR SNP was not associated with the whole-body insulinstimulated glucose disposal rate (adjusted p=0. 19). The resting metabolic rate also did not differ (p=0.94, adjusted for age, sex, fat mass and fat-free mass) between those with the A allele and those with the G allele.
Association of GCK SNPs with 24 h energy expenditure In the respiratory chamber study, participants with the A allele for the 3′UTR SNP chr7:44184184-G/A had a lower 24 h energy expenditure (by 520 kJ/day) than those with the G allele (as the difference in least square means); this was more evident in the postprandial state, as shown by the energy expenditure trajectory over the day (Fig. 2a).
Three components of 24 h energy expenditure, including the sleeping metabolic rate (60-70% of total energy expenditure), energy cost of SPA (20-30%) and thermic effect of food (awake and fed thermogenesis, 10%) were described previously [15,18]. The difference in 24 h energy expenditure between genotypes of the 3′UTR SNP is attributable to the difference in the thermic effect of food (Fig. 2b, β=520 kJ/day, p=3.39×10 −6 for 24 h energy expenditure, adjusted for age, sex, fat mass, fatfree mass and SPA), but not to SPA (p=0.16, adjusted for age and sex) or sleeping metabolic rate (adjusted for age, sex, fat mass and fat-free mass, p=0.27). To confirm this observation, we further analysed the effects of genotype on the energy expenditure trajectory during the postprandial state (daytime) vs fasting state (night-time) using a mixed model analysis. After accounting for age, sex, fat mass, fat-free mass and SPA, individuals with the A allele had a lower rate of energy expenditure during the day (β=−0.46 kJ/min, p=0.0001, from 8:00 hours on one day to 01:00 hours on the next day), compared with those with the G allele, but no difference was observed during the night (β=−0.07 kJ/min, p=0.50, from 01:00 hours to 05:00 hours). This result demonstrates that the difference in 24 h energy expenditure between genotypes is driven by the thermic effect of food in the postprandial state.
In the chamber, individuals with the A allele also had a decreased rate of carbohydrate oxidation (by 311 kJ/day) compared with those with the G allele (Fig. 2c, p=0.03, adjusted for age, sex, percentage body fat and energy balance). This was comparable with the difference in the rate of carbohydrate oxidation during a hyperinsulinaemic-euglycaemic clamp, in which a ∼70 kg (EMBS) individual with the A allele was estimated to have a lower rate of insulin-stimulated carbohydrate oxidation (by ∼406 kJ/day). However, neither the rate of lipid oxidation nor protein oxidation differed between genotypes (p=0.61 and p=0.53, respectively, adjusted for age, sex, percentage body fat and energy balance), indicating that the 3′UTR SNP affected energy expenditure primarily via carbohydrate oxidation. All eight tag SNPs were analysed for associations with metabolic traits in 415 full-heritage non-diabetic individuals (data not shown). Only the 3′UTR SNP was associated with the rate of carbohydrate oxidation and 24 h energy expenditure (Figs 1 and 2, Table 2).
Association of GCK SNPs with type 2 diabetes and BMI The eight tag SNPs were genotyped in a population-based sample of 3,604 full-heritage Pima Indians and a replication sample of 4,063 mixed-heritage American Indians for association analyses of type 2 diabetes and BMI. The 3′UTR SNP chr7:44184184-G/A had a nominal association with type 2 diabetes in full-heritage Pima Indians (Table 3; OR 1.37, 95% CI 1.06, 1.78, p=0.015, adjusted for age, sex, birth year and heritage). This SNP also had a borderline association with SNP rs4607517, which was reproducibly associated with fasting glucose concentrations in other populations [5][6][7], also had strong associations with fasting glucose concentrations in a combined analysis of 5,429 American Indians (β = 0.06 mmol/l per risk allele, p=8.7×10 −7 , adjusted for age, sex, birth year and heritage) ( Table 3). However, this SNP was not associated with type 2 diabetes (adjusted p=0.16, n= 7,667). In addition, the synonymous SNP Gly193Gly and rs193226243 were also associated with fasting glucose concentrations (β = 0.25 mmol/l, adjusted p =4.3×10 −5 ; β = 0.08 mmol/l, adjusted p=0.0003, respectively).
Beta cell function was assessed by AIR in 298 full-heritage Pima Indians with NGT. SNP rs4607517, associated with fasting glucose concentrations, was nominally associated with AIR (β=0.89 [log 10 scale], p=0.07, adjusted for age, sex, percentage body fat and rate of glucose disappearance during insulin stimulation). The promoter SNP rs193226243 was also associated with AIR (β=0.92 [log 10 scale], adjusted p= 0.02). However, the 3′UTR SNP was not associated with AIR (adjusted p=0.46).
Association of the eight tag SNPs with maximum BMI were also analysed ( Table 4). SNP rs193226243 had a nominal association with BMI in a combined analyses of 6,797 American Indians (β=0.02 (log e scale), p=0.008, adjusted for age, sex, birth year and heritage). None of other SNPs was consistently associated with BMI.
The predictive effect of carbohydrate oxidation rate on development of type 2 diabetes We further evaluated the relationship between the rate of insulin-stimulated carbohydrate oxidation and the risk of developing type 2 diabetes in 287 fullheritage Pima Indians. Of these 287 individuals with NGT who had measures of carbohydrate oxidation rate during insulin stimulation at baseline (baseline age 26.4±6.0 years) and also had follow-up data for development of diabetes (follow-up time 7.8±8.2 years), 99 (34%) developed type 2 diabetes. Figure 3 shows the Kaplan-Meier survival curve for time to type 2 diabetes onset with participants categorised as those with a higher or lower rate of insulin-stimulated carbohydrate oxidation than the mean of 1.63 mg (kg EMBS) −  . In a Cox proportional hazards analysis, a lower than the mean rate of insulin-stimulated non-oxidative glucose disposal at baseline also predicted a higher risk for developing type 2 diabetes than a rate above the mean (HRR 2.5, 95% CI 1.4, 4.5, p=0.002, adjusted for age, sex, percentage body fat, AIR and glucose oxidation rate; data not shown).
Allelic specific GCK expression To investigate whether the alleles of the 3′UTR SNP chr7:44184184-G/A differentially influence gene expression, allelic specific expression of GCK was assessed in adipose tissue biopsies from individuals heterozygous for this SNP. No difference in the allelic expression of GCK was observed in this tissue (G/A=0.993 vs expected ratio of 1, p=0.2, data not shown). Since liver and pancreas  tissue biopsies from Pima Indians were not available for study, the allelic specific expression of GCK in these tissues is not known.

Discussion
GCK is the main glucose-phosphorylating enzyme in the liver and pancreatic beta cells. It converts glucose to G6P as a first and rate-limiting step in glycolysis, which plays a part in the process of glucose oxidation [1][2][3]. Our study indicates that a novel variant in the 3′UTR of GCK, with a risk allele frequency of 0.05, is associated with a lower rate of glucose oxidation post-absorptively, during insulin stimulation and after a diet of mixed consumption, which is in agreement with the role of GCK in glycolysis. This variant was not associated with nonoxidative glucose disposal, suggesting that glucose storage (glycogen synthesis) was not affected. It is known that rare mutations in GCK occurring in MODY result from a reduced glucose-stimulated insulin secretion. Although rs193226243 and rs4607517 had borderline associations with AIR, none of the other common variants were associated with AIR. Therefore, our data suggest that common variation in GCK predominantly influences glycolysis and the rate of glucose oxidation in hepatocytes. These data are consistent with the observations that overexpression of GCK in mouse liver or rat isolated hepatocytes enhances glucose oxidation [19,20]. Nevertheless, a subtle effect of common GCK variants on beta cell function cannot be ruled out. Since the GCK variants affect the threshold for glucose sensing, the effect size and/or sample size may be too small to render a statistical difference in insulin secretion. It is also possible that the 25 g intravenous glucose bolus used in the AIR measurement may be above the threshold at which GCK exerts its effect, thus limiting a positive detection. Hepatic GCK serves as a major component of the hepatic glucose-sensing system involved in glucose synthesis, breakdown and storage. While glycolysis and glycogen synthesis pathways are activated during the postprandial state, gluconeogenesis and glycogen breakdown are involved in hepatic glucose production in the post-absorptive state. Our data indicate that the association of the 3′UTR GCK variant with the rate of basal hepatic glucose production is likely due to an association with the rate of basal glucose oxidation. However, this variant was not associated with hepatic glucose production during insulin stimulation, despite previous findings demonstrating that a GCK variant was associated with hepatic insulin resistance [21].
In addition to its pivotal role in glucose metabolism, a new role for hepatic GCK in energy metabolism has emerged in recent studies. Tsukita and co-workers reported that upregulation of hepatic GCK by high-fat diet feeding in mice suppresses brown adipose tissue (BAT) thermogenesis via leptinmediated neural signals and downregulation of uncoupling protein-1 [22,23]. This GCK-mediated liver-to-BAT neuronal relay system provides a novel mechanism in modulating obesity predisposition in mice. In this study, we report that the novel 3′UTR variant in GCK had a significant effect on 24 h energy expenditure through a change in the thermic effect of food. This effect resulted from a change in the rate of carbohydrate oxidation rather than from any apparent effect on BAT thermogenesis.
Measures of energy expenditure, thermic effect of food and substrate oxidation are predictors of weight change. In Pima Indians, individuals with lower than expected energy expenditure are at higher risk for future long-term increases in weight and fat mass [16]. The thermic effect of food is also reduced in obese individuals and predicts their future weight gain [18]. A higher rate of insulin-stimulated carbohydrate oxidation during a hyperinsulinaemic-euglycaemic clamp predicts a future weight again [24]. The rate of carbohydrate oxidation in a respiratory chamber also predicts short-term changes in body weight [14], but not long-term changes [16]. In the present study, the 3′UTR SNP, associated with rates of energy expenditure and carbohydrate oxidation, was not associated with BMI. Nevertheless, this 3′UTR SNP was associated with risk of type 2 diabetes in Pima Indians. This most likely results from the effect on carbohydrate oxidation since we found that a lower carbohydrate oxidation rate during insulin stimulation was associated with a higher risk of type 2 diabetes, independent of age, sex, percentage body fat, AIR and non-oxidative glucose disposal rate. Thus, while rare mutations in GCK cause MODY and neonatal diabetes, our data indicate that common variation in GCK with a modest effect on the rate of carbohydrate oxidation contributes to risk of type 2 diabetes. In summary, our study in individuals who had been extensively characterised for metabolic traits provides cohesive evidence to support a hepatic effect of a novel 3′UTR variant in GCK on influencing carbohydrate oxidation, energy expenditure and type 2 diabetes risk; this is consistent with the role of GCK in hepatic glycolysis and energy metabolism. However, our functional analysis of this 3′UTR SNP in adipose tissue did not support a role in allelic imbalance of GCK expression in this particular tissue. Interpretation of this negative result is unclear since this SNP could potentially affect transcriptional regulation or mRNA stability in a tissuespecific manner, and we do not have access to liver or pancreatic beta cells from Pima individuals. Alternatively, this 3′ UTR SNP might alter GCK translation via an effect on microRNA binding, or perhaps this SNP is in LD with an undiscovered functional variant. Future studies in liver or pancreatic biopsy tissues would clarify some of these possible mechanisms.