Susceptibility to type 2 diabetes may be modulated by haplotypes in G6PC2, a target of positive selection
- Nasser M. Al-Daghri†1, 2,
- Chiara Pontremoli†3,
- Rachele Cagliani3,
- Diego Forni3,
- Majed S. Alokail1, 2,
- Omar S. Al-Attas1, 2,
- Shaun Sabico1, 2,
- Stefania Riva3,
- Mario Clerici†4, 5Email author and
- Manuela Sironi†3
© The Author(s). 2017
Received: 2 November 2016
Accepted: 26 January 2017
Published: 7 February 2017
The endoplasmic reticulum enzyme glucose-6-phosphatase catalyzes the common terminal reaction in the gluconeogenic/glycogenolytic pathways and plays a central role in glucose homeostasis. In most mammals, different G6PC subunits are encoded by three paralogous genes (G6PC, G6PC2, and G6PC3). Mutations in G6PC and G6PC3 are responsible for human mendelian diseases, whereas variants in G6PC2 are associated with fasting glucose (FG) levels.
We analyzed the evolutionary history of G6Pase genes. Results indicated that the three paralogs originated during early vertebrate evolution and that negative selection was the major force shaping diversity at these genes in mammals. Nonetheless, site-wise estimation of evolutionary rates at corresponding sites revealed weak correlations, suggesting that mammalian G6Pases have evolved different structural features over time. We also detected pervasive positive selection at mammalian G6PC2. Most selected residues localize in the C-terminal protein region, where several human variants associated with FG levels also map. This region was re-sequenced in ~560 subjects from Saudi Arabia, 185 of whom suffering from type 2 diabetes (T2D). The frequency of rare missense and nonsense variants was not significantly different in T2D and controls. Association analysis with two common missense variants (V219L and S342C) revealed a weak but significant association for both SNPs when analyses were conditioned on rs560887, previously identified in a GWAS for FG. Two haplotypes were significantly associated with T2D with an opposite effect direction.
We detected pervasive positive selection at mammalian G6PC2 genes and we suggest that distinct haplotypes at the G6PC2 locus modulate susceptibility to T2D.
KeywordsG6PC2 G6PC G6PC3 Natural selection Association analysis Type 2 diabetes
The endoplasmic reticulum enzyme glucose-6-phosphatase catalyzes the hydrolysis of glucose-6-phosphate (G6P) to glucose and inorganic phosphate. The enzyme is part of a multicomponent integral membrane system that includes the catalytic subunit (G6PC, hereafter referred to as G6Pase) as well as transporters for glucose-6-phosphate, inorganic phosphate, and glucose [1, 2]. G6Pase catalyzes the common terminal reaction in the gluconeogenic and glycogenolytic pathways, resulting in the release of glucose into the bloodstream . These results led to the identification of G6Pase as a key player in glucose homeostasis.
In most mammals, different G6PC subunits are encoded by three paralogous genes (G6PC, G6PC2, and G6PC3), usually referred to as the G6PC gene family [1, 2]. The protein products of the three genes display moderate sequence identity and a common topological organization with nine transmembrane domains and intralumenal catalytic residues .
G6PC is mainly expressed in the liver and kidney and at lower levels in the intestine and pancreatic islets, and has a critical function in maintaining euglycemia in fasting conditions [1, 2]. In humans, mutations in the gene cause glycogen storage disease type Ia (GSD1A), which results in severe hypoglycemia and glycogen accumulation-associated hepatomegaly, as well as growth retardation, lactic acidemia, hyperlipidemia, hyperuricemia, and increased incidence of hepatic adenomas [1, 2]. Mutations in G6PC3 are also associated with pathology in humans. Thus, although the gene is ubiquitously expressed, its function is particularly important in white blood cells, and G6PC3 deficiency causes autosomal recessive severe congenital neutropenia type 4 (SCN4) [1, 2]. SCN4 patients are particularly susceptible to bacterial infections and may display additional non immunologic symptoms. Conversely, in both humans and in the knock-out mouse model, G6PC3 only marginally contributes to the regulation of blood glucose levels or hepatic glycogen content [1, 2]. Finally, G6PC2 is specifically expressed in pancreatic islets where its function is still incompletely understood [1, 2]. g6pc2 −/− mice display a reduction in blood glucose levels after a 6 h fast, whereas plasma insulin and glucagon concentrations are unaffected [1, 2]. These data led to the hypothesis that G6PC2 regulates the glycolytic flux by hydrolyzing G6P, thereby opposing the action of glucokinase. G6PC2 and glucokinase are, therefore suggested to function as beta islet glucose sensors [1, 2]. In humans, common and rare variants in G6PC2 have been associated with fasting glucose (FG) levels and with decreased insulin secretion during glucose tolerance tests [3–9]. This observation led to the suggestion that G6PC2 may also regulate the pulsatility of insulin secretion [1, 2].
Variation in FG is clinically important in humans, as it is associated with the risk of developing type 2 diabetes (T2D) and ischemic heart disease [10, 11] as well as being an important determinant of offspring birth weight in pregnant women .
In humans and other mammals, FG levels are influenced by the feeding status. Prolonged fasting causes a reduction in blood glucose levels, which can result in life-threatening hypoglycemia; the gluconeogenic pathway is the major contributor to the maintenance of glucose levels during fasting and starvation . Mammals display a wide variety of diets, different lifestyles (that may or may not include recurrent prolonged fasts), and distinct energy requirements. These characteristics influence the ability of a species to resist prolonged fasting , a situation that is common in nature and that is likely to exert a strong selective pressure. It is thus conceivable that genes involved in the regulation of FG have been targeted by positive (or diversifying) selection during mammalian evolution. Indeed, positive selection was previously demonstrated to act on genes with a role in carbohydrate absorption and digestion in mammals [14, 15]. In humans, aside from the textbook example of lactase persistence , signals of diet-driven selection include variants in genes involved in starch and sucrose metabolism [15, 17], copy number variation at genes encoding salivary amylase (AMY1) , as well as polymorphisms in genes that may be associated with the consumption of cooked food . In fact, humans likely underwent several dietary shifts associated with cultural innovations such as the use of fire for cooking (likely predating the split of modern humans from Neanderthals/Denisovans) [19, 20], the exploitation of starch-rich plant underground storage organs , and the agricultural revolution. Because these cultural changes modified diet composition and caloric intake, genes involved in glucose homeostasis, such as G6PC genes, represent likely target of positive selection in humans.
Herein we use both inter- and intra-species comparisons to analyze the evolution of the three G6Pase genes in mammals and human populations. We also perform an association study to assess the role of G6PC2 variants in T2D susceptibility in a population with high incidence of metabolic disorders.
Evolutionary origin of the G6PC gene family
We first investigated the evolutionary origin of the three mammalian G6PC paralogs. Analysis of a gene gain/loss tree of 70 animal species through the Ensembl Compara utility [22, 23] indicated that a single G6PC gene is present in the Drosophila genome, whereas lamprey (Petromyzon marinus, Cyclostomata) displays two genes and most bony fishes, birds, reptiles, amphibians and mammals have at least three paralogs. Possibly due to gene loss, no G6PC gene is described in the two Tunicata genomes included in the Ensembl Compara dataset.
Analysis of the G6PC phylogeny indicated that an initial duplication event in the lineage basal to all vertebrates originated G6PC3 and the G6PC/G6PC2 ancestor. In lamprey, one of the two G6PC sequences clusters with G6PC3 proteins, whereas the other is basal to G6PC2 and G6PC (Fig. 1), suggesting that the duplication events that originated G6PC and G6PC2 occurred after the split of gnathostomes and cyclostomes but before the divergence of cartilaginous and bony fishes, as the three Callorhinchus milii sequences (the elephant shark) indicate (Fig. 1).
Evolutionary analysis of the glucose-6-phosphatase (G6PC) catalytic subunit gene family in mammals
Average non-synonymous/synonymous substitution rate ratio (dN/dS) and percentage of negatively selected sites fot the three G6Pase genes
Protein size (amino acids)
N° of species
Average dN/dS (95% confidence intervals)
% of FEL negatively selected sites
0.167 (0.156, 0.178)
0.206 (0.191, 0.222)
0.159 (0.147, 0.171)
A common expectation is that mutations at highly constrained codons are more likely to disrupt protein structure/function and, therefore, to cause disease. To date, 57 independent GSD1A missense mutations involving 47 unique codons have been reported. We observed that codons that carry at least one missense mutation are significantly more likely to show statistical evidence of negative selection (FEL p value < 0.1) than codons where no mutation has been described (Fisher’s Exact Test, two tailed, p = 0.044, odds ratio = 2.19, 95% confidence intervals: 0.99–5.34). The same calculation was not performed for G6PC3 mutations, as too few of such mutations are actually known (number of mutated codons = 9, seven of which negatively selected).
Positive selection at the mammalian G6PC2 gene
Likelihood ratio test statistics for models of variable selective pressure among sites in G6PC2
Codon frequency model
Degrees of freedom
% of sites (average dN/dS)
Positively selected sites
M1a vs M2aa
M7 vs M8b
G137 (BEB, REL, FEL), A297 (BEB), L298 (BEB, REL, FEL), E316 (BEB), G351 (BEB, REL)
M8ac vs M8
M1a vs M2a
M7 vs M8
M8a vs M8
Evolutionary analysis of G6Pase genes in humans and great apes
We next applied a population genetics-phylogenetics approach to study the evolution of G6Pase genes in the human, chimpanzee, and gorilla lineages. Specifically, we ran the gammaMap program  that jointly uses intra-specific variation and inter-specific diversity to estimate the distribution of fitness effects (i.e. population-scaled selection coefficients, γ) along coding regions. gammaMap categorizes codons into 12 classes of γ, ranging from strongly beneficial (γ = 100) to inviable (γ = −500); a γ equal to 0 indicates neutrality. The overall distribution of selection coefficients indicated that G6PC evolved under strong purifying selection in all lineages (median γ < −10, Fig. 2b). This was also the case for G6PC2 in non-human primates (median γ = −100), whereas the human gene showed weaker constraint (Fig. 2b). Finally, the distribution of fitness effects for G6PC3 was very different in distinct lineages. In fact, the codon distribution was almost homogeneous across the range of γ values in humans and chimpanzees, although the median remained below 0. In contrast, the gorilla lineage showed evidence of strong purifying selection (Fig. 2b). We thus assessed whether this pattern may derive from a relaxation of constraint in humans and chimpanzees. To test this possibility we applied the RELAX methodology  to the G6PC3 primate phylogeny (Fig. 2c). Results were consistent with relaxed selection on the human/chimpanzee branches (p = 0.037, k = 0), but not on the gorilla lineage (p = 0.958, k = 1.05) (Fig. 2c). The same analysis for the human G6PC2 branch revealed no relaxation of selective pressure (p = 0.866, k = 1.21). gammaMap also identified two positively selected codons (cumulative probability > 0.80 of γ ≥ 1) for human G6PC2 (Fig. 2, Additional file 1: Table S4). One selected codon was also identified for human G6PC3 (site 243), whereas no positively selected sites were detected for G6PC in any lineage.
Evolutionary analysis in human populations
We finally investigated whether positive selection acted on G6Pase genes during the recent evolutionary history of human populations. Using the 1000 Genomes Phase 1 data for Yoruba, European, and Chinese we calculated pairwise FST , an estimate of population genetic differentiation. We also performed the DIND (Derived Intra-allelic Nucleotide Diversity) and iHS (integrated haplotype score) tests [35, 36] for all SNPs mapping to these genes. Statistical significance (in terms of percentile rank) for the FST statistic and for the DIND test was obtained by deriving empirical distributions. For the iHS test, absolute values higher than 2 were considered as significant . No SNP in any G6Pase gene reached statistical significance (rank > 0.95) for both FST and for the DIND tests, and none had an |iHS| higher than 2. Overall, these results indicate that no variant/haplotype can be confidently called as positively selected. Likewise, nucleotide diversity (calculated as θW and π [37, 38]) for the entire gene regions was unexceptional if compared to those calculated for a reference set of 2000 genes. We conclude that G6Pase genes did not represent selection targets in recent human history.
Association of G6PC2 variants with T2D
Association of G6PC2 variants with T2D
Minor allele freq (%)
Corrected p value
OR (IC 95%)
Corrected p value
OR (IC 95%)
Conditional on rs560887
rs560887, intronic, (GWAS)
G6PC2 haplotype analysis
Frequency in T2D (%)
Frequncy in unaffected (%)
Association p value
rs560887 | rs492594| rs2232328
Finally, to assess the effect of rare and common G6PC2 variants on T2D risk, we applied a SNP-set based method, the Sequence Kernel Association Test (SKAT) . SKAT was run either by inclusion of all variants identified through re-sequencing (n = 13, Fig. 3, Additional file 1: Tables S6 and S7) or by limiting analysis to missense SNPs plus the GWAS variant (rs560887). No significant association was detected. However, as for single-variant associations, the power of SKAT is limited when small samples are analyzed .
In this study we have analyzed the evolutionary history of three genes (G6PC, G6PC2 and G6PC3) encoding the catalytic subunits of glucose-6-phosphatase, a central enzyme for glucose homeostasis. The analysis was motivated by the well-accepted concept that the availability of food resources is a driver of pivotal importance in the evolution in mammals and that, in natural settings, most mammals commonly face prolonged fasting and/or starvation . Consequently, homeostatic mechanisms that sense plasma glucose levels and modulate them in response to the feeding status are expected to represent natural selection targets.
Commonly, positive and negative selection act in concert on the same protein-coding gene. In fact, due to structural and functional constraints, most amino acid replacements are deleterious and are eliminated by negative selection. Conversely, at a minority of sites, amino acid replacements may be favored because, without impairing protein function, they confer new advantageous properties . In line with this view, we found all G6Pase genes to display an overall dN/dS lower than 1, indicating a preponderance of negative selection. Recent evidence showed that structural and folding requirements (i.e. the ability of a protein to fold properly and stably) represent major determinants of the evolutionary rate at protein sites . The 3D structures of mammalian G6Pases has not been solved and we could not therefore assess whether among-site variation in evolutionary rates is correlated with parameters such as solvent accessibility or packing density . Nonetheless, we reasoned that because the three proteins share considerable identity in terms of amino acid sequence and the same topological organization , they should also display a similar 3D structure and, consequently, corresponding residues should display similar evolutionary rates. In fact, this was only partially true, as the correlation of dN-dS at corresponding sites were either weak or non-significant. This suggests that, despite a similar membrane topology ad the maintenance of the catalytic function, mammalian G6Pases have evolved different structural features over time. Indeed, the three genes have been diverging for a long time, as the duplications that originated the three mammalian paralogs occurred during early vertebrate evolution. It is generally accepted  that two whole genome duplication events occurred in the lineage basal to all vertebrates, before the divergence of gnathostomes and cyclostomes, although some authors favored a model with a single whole genome duplication . It is thus possible that G6PC3 and the G6PC/G6PC2 ancestor originated and were fixed after whole genome duplication(s) in the ancestral vertebrate. However, the basal position of one lamprey sequence with respect to gnathostome G6PC and G6PC2 proteins suggests that the duplication event that originated the two genes occurred after the gnathostome/cyclostome split. After gene duplications, gene losses occurred in several species or lineages; for instance most marsupials and the platypus only have one G6PC gene. Additional G6PC duplications also occurred during vertebrate evolution; several bony fishes have 4 G6PC paralogs, possibly as a results of a whole genome duplication that occurred in the ancestor of teleosts . A similar observation was reported for the rainbow trout, a glucose-intolerant fish, which displays 5 G6PC genes possibly fixed in this species after the salmonid-specific whole genome duplication . Overall, these observations indicate that the G6PC gene family is highly dynamic and gene maintenance or loss in some lineages may be related to specific feeding needs or strategies.
In line with this view, we detected pervasive positive selection at mammalian G6PC2 genes. Most residues targeted by selection are located in the C-terminal protein region, which is also subject to strong negative selection. Because of the role of G6PC2 as a glucose sensor, it is possible to speculate that adaptive changes in distinct mammals relate to trophic strategies including diet, hybernation, and feeding behavior. Interestingly, positively selected sites in the human G6PC2 gene were detected as well. It is worth mentioning that the two selected residues are fixed or almost fixed in human populations; checking against the genome sequences of archaic hominins indicated that the C46 and A119 variant were already present in the genomes of Neandertals and Denisovans [50, 51]. These observations suggest that, as for other variants in metabolic genes , these changes were not driven to high frequency in humans as an adaptation to the dietary shift determined by agriculture. Indeed, population genetics analysis of modern human populations detected no recent selective event.
Unexpectedly, given its association with a human disease, two different analyses indicated that G6PC3 genes have experienced a relaxation of selective pressure in the human and chimpanzee lineages. We note, however, that this finding does not imply that relaxed constraints are observed at all sites in the protein. Conversely, in humans this effect is driven by 4 nonsynonymous substitutions (either fixed or polymorphic relative to the common ancestor of Hominidae), including the positively selected 243 site, in the absence of synonymous substitution. Three of these changes are clustered in ~60 amino acid region (residues 216–275) suggesting that, for unknown reasons, this protein portion is tolerant to change in humans. To date, no SNC4 missense mutation has been described at these sites.
Among the three G6Pase genes, mutations in G6PC2 have never been associated with a Mendelian human disease. This is in line with the mild phenotype of the KO mouse model, as well as with the observation that G6PC2 is a pseudogene in rats. Moreover, recent functional data indicated that coding variants that reduce the expression of G6PC2, most likely by impairing its folding, segregate at appreciable frequency in human populations . Notably, variants in G6PC2 have been consistently associated with FG levels, whereas their contribution to T2D risk remains controversial. In particular, the rs560887 SNP is one of the strongest signals associated to FG (and related traits), and one of the most commonly replicated in large-scale analyses [3–6, 52–54]. Moreover, the variant was shown to be functional and to modulate G6PC2 pre-mRNA splicing . Although this latter finding does not necessarily imply that rs560887 is the causal variant, the effect of the glucose-raising allele (C) on increased splicing efficiency is suggestive . However, distinct studies found either no association of rs560887 with T2D risk  or indicated a weak protective effect of the glucose-increasing allele [6, 41]. Recently, Mahajan and coworkers reported a glucose-increasing effect of the common V219 (rs492594-G) allele that modestly increases the risk of T2D as well . The authors suggested that association analysis for G6PC2 should be performed through haplotype reconstruction as multiple rare and common variants independently affect FG levels, and the direction of effect for rs492594 is reversed when analysis is conditioned on rs560887 . Nonetheless, most large-scale analyses of T2D susceptibility performed single variant association tests, rather than haplotype inference, leaving the role of G6PC2 in T2D partially unexplored.
Our sequencing analysis in the Saudi sample was motivated by the high prevalence of T2D in this population. The frequency of rare variants was not different in T2D and HC, but the small sample size is not well suited to this type of analysis. Haplotype analysis with common variants detected two haplotypes that associated with T2D susceptibility in Saudi subjects. The haplotypes include the rs2232328 (S342C) variant, that is not covered in exome chip arrays and was thus not analyzed in recent association studies of G6PC2 variants for FG levels [39, 40]. In a genome-wide meta-analysis , rs2232328 showed a strong association with FG (p value adjusted for BMI = 5.1 × 10−16), which is likely independent of the lead variant rs560887, as their LD is low in all populations (r2 < 0.05) (http://analysistools.nci.nih.gov/LDlink/). The functional effect of the S342C substitution is presently unknown. Codon 342 is negatively selected in mammals and located in a highly constrained region; indeed, a cystein residue was present in all analyzed mammals with the only exception of macaques (Additional file 1: Figure S1). These observations suggest that the derived S342 allele impairs G6PC2 function. Surprisingly, though, the V219 allele which also involves a negatively selected site and represents the ancestral state conserved in all mammals (with the only exception of the tree shrew), was recently shown to result in reduced function . Indeed, G6PC2 molecules carrying the V219 allele are expressed at lower abundance due to proteasomal degradation . This observation indicates that the functional effect of G6PC2 variants is difficult to predict, and in the case of the S342 substitution will need experimental testing.
The data we report herein, although preliminary, may help reconcile the contrasting results obtained for rs560887 on T2D risk, as its effect might depend on haplotype context and may vary in different populations depending on LD between rs560887 and other functional variants.
Clearly, further studies will be necessary to confirm the role of G6PC2 variants on T2D susceptibility. First, the size of the Saudi sample is small and the associations we detected are weak, thus requiring validation in an independent larger sample. Second, variants in the 5′ region of G6PC2 (rs13387347, rs1402837) and in the intergenic spacer downstream the transcription end site of the gene (rs563694) were also associated with FG [4, 55, 56]. These variants possibly contribute independently to FG levels and show variable levels of LD with the SNPs we analyzed. Because the focus of our work was on coding missense variants, we did not analyze these SNPs. However, they may contribute to T2D susceptibility either alone or in combinations with coding variants, warranting their inclusion in future efforts aimed at assessing the contribution of G6PC2 genetic variability to T2D risk.
In conclusion, we detected pervasive positive selection at mammalian G6PC2 genes, with almost all selected sites located in the C-terminal portion of the protein.
We then investigated a possible role for G6PC2 variants in modulating the susceptibility to T2D in subjects from Saudi Arabia. We detected two haplotypes, one predisposing and one protective, significantly associated with T2D. These preliminary results suggest that distinct G6PC2 haplotypes modulate susceptibility to T2D.
Phylogenetic analysis in metazoans
Protein sequences of G6PC genes for 65 animal species were retrieved from the Ensembl Compara database (Additional file 1: Table S1). The genomes of the following metazoans were searched for G6PC orthologs and paralogs: Strongylocentrotus purpuratus, Aplysia californica, Callorhinchus milii, Saccoglossus kowalevskii, Limulus polyphemus. Searches were performed using BLASTp using the three human G6PC proteins as queries, as well as the two lamprey proteins and the single protein of sea urchin. All hits corresponded to predicted proteins derived from genomic sequences.
The genomes of three Cephalochordata (Branchiostoma lanceolatum, Branchiostoma belcheri, and Asymmetron lucayanum) was also searched for the presence of G6PC genes but no hit was obtained.
A maximum likelihood phylogenetic tree of 188 G6PC proteins was constructed using RAxML v8.2.9  with 100 bootstrap replicates and the best protein substitution model automatically determinated by the software.
Evolutionary analysis in mammals
Available mammalian sequences for G6PC, G6PC2 and G6PC3 were retrieved from the NCBI database (http://www.ncbi.nlm.nih.gov/). A list of species is available as Additional file 1: Table S2. DNA alignments were performed with the RevTrans 2.0 utility (http://www.cbs.dtu.dk/services/RevTrans/, MAFFT v6.240 as an aligner) , which uses the protein sequence alignment as a scaffold for constructing the corresponding DNA multiple alignment. All alignments were screened for the presence of recombination using GARD (Genetic Algorithm Recombination Detection) , a Genetic Algorithm implemented in the HYPHY suite . Gene trees were generated by maximum-likelihood using phyML with a maximum-likelihood approach, a General Time Reversible (GTR) model plus gamma-distributed rates and 4 substitution rate categories .
The SLAC (Single-Likelihood Ancestor Counting) and FEL (Fixed Effects Likelihood) methods from the HYPHY package were used to calculate the overall dN/dS, to identify negatively selected sites (FEL significance cut-off = 0.1) and for calculating dN-dS (rate of nonsynonymous changes-rate of synonymous changes) at each site .
The site models implemented in PAML were developed to detect positive selection affecting only a few aminoacid residues in a protein. To detect selection, site models that allow (M2a, M8) or disallow (M1a, M7 and M8a) a class of sites to evolve with ω > 1 were fitted to the data using the F3x4 (codon frequencies estimated from the nucleotide frequencies in the data at each codon site) and the F61 (frequencies of each of the 61 non-stop codons estimated from the data) codon frequency model. Positively selected sites were identified using the Bayes Empirical Bayes (BEB) analysis (with a cut-off of 0.95). BEB calculates the posterior probability that each codon is from the site class of positive selection (under model M8) . The REL (Random Effects Likelihood)  and FEL (with the default cutoff of 0.1) tools were also applied to identify positively selected sites. REL models variation in nonsynonymous and synonymous rates across sites according to a predefined distribution, with the selection pressure at an individual site inferred using an empirical Bayes approach; FEL directly estimates nonsynonymous and synonymous substitution rates at each site .
Tests for potential-relaxed selection of G6PC2 and G6PC3 genes in primates were performed using the hypothesis testing framework in RELAX from the HYPHY package . RELAX calculates a selection intensity parameter, k, by taking into account that relaxation will exert different effects on sites subjected to purifying selection (ω < 1) and sites subjected to positive selection (ω > 1). Relaxation will move ω toward 1 for both categories. RELAX tests whether selection is relaxed or intensified on a subset of test branches compared with a subset of reference branches in a predefined tree. In the null model, the selection intensity is constrained to 1 for all branches, whereas in the alternative model k is allowed to differ between reference and test groups. The selection on test branches is intensified or relaxed compared with background branches when k > 1 or k < 1, respectively.
Positive selection in Homininae
For gammaMap  analysis, genotype data from the phase 1 of the 1000 Genomes Project were retrieved from the dedicated website ; we retrieved SNP information for the three human populations: African (Yoruba), European, and Chinese. For the chimpanzee and gorilla analyses, genotype information were retrieved from  for 25 and 27 individuals, respectively.
Ancestral sequences were reconstructed by parsimony from the human, chimpanzee, orangutan and macaque sequences.
Analysis was performed assuming θ (neutral mutation rate per site), k (transitions/transversions ratio), and T (branch length) to vary among genes following log-normal distributions. For p (the probability that adjacent codons share the same population-scaled selection coefficient) we assumed a uniform distribution. For each gene we set the neutral frequencies of non-STOP codons (1/61). For selection coefficients, we considered a uniform Dirichlet distribution with the same prior weight for each selection class. For each gene we run 500,000 iterations with a 20,000 iteration burn-in and a thinning interval of 10 iterations.
G6PC and G6PC3 mutations
The list of G6PC and G6PC3 mutations was obtained from the Human Gene Mutation Database (HGMD, http://www.hgmd.cf.ac.uk/ac/) and the ClinVar database (http://www.ncbi.nlm.nih.gov/clinvar/). Only missense mutations were included in the analyses.
Population genetics analyses
Genotype information was obtained for the 3 genes (G6PC, G6PC2 and G6PC3); in particular, three human populations with different ancestry were analyzed: Europeans (CEU), Africans (Yoruba,YRI), and East Asians (Han Chinese in Bejing, CHB). A control set of ~2,000 randomly selected genes was used as a reference set (hereafter referred to as control set). These gene were selected to be longer than 5000 bp and have more than 80% human-outgroup (chimpanzee, orangutan or macaque genomes) aligning bases; orthologous regions in the outgroups were retrieved using the LiftOver tool. These data were used to calculate θW , π  and their relative distruibutions.
The pairwise FST  and the DIND (Derived Intra-allelic Nucleotide Diversity)  test were calculated for all SNPs mapping to the analyzed genes, as well as for SNPs mapping to the control set. FST values are not independent from allele frequencies, so we binned variants in 50 classes based on the minor allele frequency (MAF) and calculated F ST empirical distribution for each MAF class using the control set data. The same procedure was applied for the DIND test; thus, we calculated statistical significance by obtaining an empirical distribution of DIND values for variants located within control genes; in particular, the DIND test was calculated using a constant number of 40 flanking variants (20 upstream and 20 downstream), as previously described . DIND values for the three human populations were binned in 100 derived allele frequency (DAF) classes, and for each class the distributions were calculated. As suggested , for values of iπD = 0 we set the DIND value to the maximum obtained over the corresponding class plus 20.
The iHS statistic was calculated as previously described  for all variants mapping to G6PC genes. Specifically, the iHS value was calculated using all SNPs surrounding each variant in a 5 kb region.
Human subjects, genotyping and statistical analysis
A total of 562 subjects from the Biomarker Screening in Riyadh Project (RIYADH COHORT) were enrolled (Additional file 1: Table S5). Diagnosis of T2D was based on World Health Organization proposed cut-off: fasting plasma glucose > or = 7.0 mmol/L or 126 mg/dl. Subjects with medical complications (coronary artery disease, nephropathy, and end stage renal disease or liver disease) were excluded and a similar percentage of males and females was enrolled among T2D patients and controls. Anthropometry included measurement of height (to the nearest 0.5 cm) and weight (to the nearest 0.1 kg); BMI was calculated as kg/m2. According to the WHO criteria, individuals were classified as non obese if their BMI was < 30 kg/m2. Written consent was obtained from all participants, and ethical approval was granted by the Ethics Committee of the College of Science Research Center, King Saud University, Riyadh, Kingdom of Saudi Arabia (KSA).
The two terminal exons of G6PC2 were resequenced through PCR amplification and direct sequencing (primer sequences are available upon request). PCR products were treated with ExoSAP-IT (USB Corporation Cleveland Ohio, USA), directly sequenced on both strands with a Big Dye Terminator sequencing Kit (v3.1 Thermo Fisher Scientific), and run on an Applied Biosystems ABI 3130 XL Genetic Analyzer. Sequences were assembled using DNA Baser Sequence Assembler version 4.10. A summary of all variants identified through resequencing is available (Fig. 3, Table 3, Additional file 1: Table S6 and S7).
Genetic association was investigated by logistic regression with age, sex and BMI as covariates, conditioning or not on the GWAS variant (rs560887). Analyses were performed using PLINK .
The SKAT test is implemented in the SKAT R package . The SKAT_commonRare function was used to combine the effect of common and rare variants. The suggested threshold (1/√2n, where n is the number of subjects) was used to define rare variants . Analyses were performed either by deriving variant weights from a beta density function  or by using weights based on the minor allele frequency in the analyzed populations. As in PLINK analysis, age, sex and BMI were entered as covariates.
1000 Genomes Pilot Project
Bayes Empirical Bayes
Body Mass Index
(Chinese plus Japanese)
Derived Allele Frequency
Derived Intra-allelic Nucleotide Diversity
Non-synonymous substitution/synonymous substitution rate
Fixed Effects Likelihood
Genetic algorithm recombination detection
Genome-Wide Association Study
likelihood ratio test
Minor Allele Frequency
Mixed Effects Model of Evolution
Phylogenetic Analysis by Maximum Likelihood
Random Effects Likelihood
Single-likelihood ancestor counting
Single Nucleotide Polymorphism
Type 2 Diabetes
This work is partially supported by Prince Mutaib Bin Abdullah Chair for Biomarkers of Osteoporosis of King Saud University and by the Italian Ministry of University and Research.
CP is supported by a fellowship of the Doctorate School of Molecular Medicine, University of Milan.
Availability of data and materials
Multiple protein alignment of G6PC2 genes is provided as Additional file 1: Figure S1, where positively selected sites are also highlighted.
Primer sequences are available upon request.
Alignments, trees and input files for gammaMap, as well as files used for PLINK and SKAT tests, were deposited in the Dryad repository. Dryad DOI: doi:10.5061/dryad.pk50t
CP, MS, RC, DF, and SR performed the experiments and analysed the data. MC and MS designed the study and contributed to writing the manuscript, with inputs from all authors. NM AD contributed to the study design and to writing the manuscript. NM AD, MSA, OS AA, and SS are responsible for following the patient cohorts and for collecting and cataloguing the samples. All authors read and approved the manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Written consent was obtained from all participants, and ethical approval was granted by the Ethics Committee of the College of Science Research Center, King Saud University, Riyadh, Kingdom of Saudi Arabia (KSA).
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Marcolongo P, Fulceri R, Gamberucci A, Czegle I, Banhegyi G, Benedetti A. Multiple roles of glucose-6-phosphatases in pathophysiology: state of the art and future trends. Biochim Biophys Acta. 2013;1830(3):2608–18.View ArticlePubMedGoogle Scholar
- O’Brien RM. Moving on from GWAS: functional studies on the G6PC2 gene implicated in the regulation of fasting blood glucose. Curr Diab Rep. 2013;13(6):768–77.View ArticlePubMedPubMed CentralGoogle Scholar
- Bouatia-Naji N, Rocheleau G, Van Lommel L, Lemaire K, Schuit F, Cavalcanti-Proenca C, et al. A polymorphism within the G6PC2 gene is associated with fasting plasma glucose levels. Science. 2008;320(5879):1085–8.View ArticlePubMedGoogle Scholar
- Chen WM, Erdos MR, Jackson AU, Saxena R, Sanna S, Silver KD, et al. Variations in the G6PC2/ABCB11 genomic region are associated with fasting glucose levels. J Clin Invest. 2008;118(7):2620–8.PubMedPubMed CentralGoogle Scholar
- Reiling E, van’t Riet E, Groenewoud MJ, Welschen LM, van Hove EC, Nijpels G, et al. Combined effects of single-nucleotide polymorphisms in GCK, GCKR, G6PC2 and MTNR1B on fasting plasma glucose and type 2 diabetes risk. Diabetologia. 2009;52(9):1866–70.View ArticlePubMedPubMed CentralGoogle Scholar
- Dupuis J, Langenberg C, Prokopenko I, Saxena R, Soranzo N, Jackson AU, et al. New genetic loci implicated in fasting glucose homeostasis and their impact on type 2 diabetes risk. Nat Genet. 2010;42(2):105–16.View ArticlePubMedPubMed CentralGoogle Scholar
- Baerenwald DA, Bonnefond A, Bouatia-Naji N, Flemming BP, Umunakwe OC, Oeser JK, et al. Multiple functional polymorphisms in the G6PC2 gene contribute to the association with higher fasting plasma glucose levels. Diabetologia. 2013;56(6):1306–16.View ArticlePubMedPubMed CentralGoogle Scholar
- Rose CS, Grarup N, Krarup NT, Poulsen P, Wegner L, Nielsen T, et al. A variant in the G6PC2/ABCB11 locus is associated with increased fasting plasma glucose, increased basal hepatic glucose production and increased insulin release after oral and intravenous glucose loads. Diabetologia. 2009;52(10):2122–9.View ArticlePubMedGoogle Scholar
- Heni M, Ketterer C, Hart LM, Ranta F, van Haeften TW, Eekhoff EM, et al. The impact of genetic variation in the G6PC2 gene on insulin secretion depends on glycemia. J Clin Endocrinol Metab. 2010;95(12):E479–84.View ArticlePubMedGoogle Scholar
- Tirosh A, Shai I, Tekes-Manova D, Israeli E, Pereg D, Shochat T, Kochba I, Rudich A, Israeli Diabetes Research Group. Normal fasting plasma glucose levels and type 2 diabetes in young men. N Engl J Med. 2005;353(14):1454–62.View ArticlePubMedGoogle Scholar
- Bjornholt JV, Erikssen G, Aaser E, Sandvik L, Nitter-Hauge S, Jervell J, Erikssen J, Thaulow E. Fasting blood glucose: an underestimated risk factor for cardiovascular death. Results from a 22-year follow-up of healthy nondiabetic men. Diabetes Care. 1999;22(1):45–9.View ArticlePubMedGoogle Scholar
- Breschi MC, Seghieri G, Bartolomei G, Gironi A, Baldi S, Ferrannini E. Relation of birthweight to maternal plasma glucose and insulin concentrations during normal pregnancy. Diabetologia. 1993;36(12):1315–21.View ArticlePubMedGoogle Scholar
- McCue MD. Starvation physiology: reviewing the different strategies animals use to survive a common challenge. Comp Biochem Physiol A Mol Integr Physiol. 2010;156(1):1–18.View ArticlePubMedGoogle Scholar
- Axelsson E, Ratnakumar A, Arendt ML, Maqbool K, Webster MT, Perloski M, et al. The genomic signature of dog domestication reveals adaptation to a starch-rich diet. Nature. 2013;495(7441):360–4.View ArticlePubMedGoogle Scholar
- Pontremoli C, Mozzi A, Forni D, Cagliani R, Pozzoli U, Menozzi G, Vertemara J, Bresolin N, Clerici M, Sironi M. Natural Selection at the Brush-Border: Adaptations to Carbohydrate Diets in Humans and Other Mammals. Genome Biol Evol. 2015;7(9):2569–84.View ArticlePubMedPubMed CentralGoogle Scholar
- Tishkoff SA, Reed FA, Ranciaro A, Voight BF, Babbitt CC, Silverman JS, Powell K, et al. Convergent adaptation of human lactase persistence in Africa and Europe. Nat Genet. 2007;39(1):31–40.View ArticlePubMedGoogle Scholar
- Hancock AM, Witonsky DB, Ehler E, Alkorta-Aranburu G, Beall C, Gebremedhin A, et al. Colloquium paper: human adaptations to diet, subsistence, and ecoregion are due to subtle shifts in allele frequency. Proc Natl Acad Sci U S A. 2010;107 Suppl 2:8924–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Perry GH, Dominy NJ, Claw KG, Lee AS, Fiegler H, Redon R, et al. Diet and the evolution of human amylase gene copy number variation. Nat Genet. 2007;39(10):1256–60.View ArticlePubMedPubMed CentralGoogle Scholar
- Carmody RN, Dannemann M, Briggs AW, Nickel B, Groopman EE, Wrangham RW, Kelso J. Genetic Evidence of Human Adaptation to a Cooked Diet. Genome Biol Evol. 2016;8(4):1091–103.View ArticlePubMedPubMed CentralGoogle Scholar
- Henry AG, Brooks AS, Piperno DR. Microfossils in calculus demonstrate consumption of plants and cooked foods in Neanderthal diets (Shanidar III, Iraq; Spy I and II, Belgium). Proc Natl Acad Sci U S A. 2011;108(2):486–91.View ArticlePubMedGoogle Scholar
- Laden G, Wrangham R. The rise of the hominids as an adaptive shift in fallback foods: plant underground storage organs (USOs) and australopith origins. J Hum Evol. 2005;49(4):482–98.View ArticlePubMedGoogle Scholar
- Vilella AJ, Severin J, Ureta-Vidal A, Heng L, Durbin R, Birney E. EnsemblCompara GeneTrees: Complete, duplication-aware phylogenetic trees in vertebrates. Genome Res. 2009;19(2):327–35.View ArticlePubMedPubMed CentralGoogle Scholar
- Herrero J, Muffato M, Beal K, Fitzgerald S, Gordon L, Pignatelli M, et al.: Ensembl comparative genomics resources. Database (Oxford) 2016, 2016:10.1093/database/bav096. Print 2016.
- Martin CC, Bischof LJ, Bergman B, Hornbuckle LA, Hilliker C, Frigeri C, et al. Cloning and characterization of the human and rat islet-specific glucose-6-phosphatase catalytic subunit-related protein (IGRP) genes. J Biol Chem. 2001;276(27):25197–207.View ArticlePubMedGoogle Scholar
- Kosakovsky Pond SL, Posada D, Gravenor MB, Woelk CH, Frost SD. Automated phylogenetic detection of recombination using a genetic algorithm. Mol Biol Evol. 2006;23(10):1891–901.View ArticlePubMedGoogle Scholar
- Kosakovsky Pond SL, Frost SD. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Mol Biol Evol. 2005;22(5):1208–22.View ArticlePubMedGoogle Scholar
- Sironi M, Cagliani R, Forni D, Clerici M. Evolutionary insights into host-pathogen interactions from mammalian sequence data. Nat Rev Genet. 2015;16(4):224–36.View ArticlePubMedGoogle Scholar
- Yang Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol. 2007;24(8):1586–91.View ArticlePubMedGoogle Scholar
- Yang Z. PAML: a program package for phylogenetic analysis by maximum likelihood. Comput Appl Biosci. 1997;13(5):555–6.PubMedGoogle Scholar
- Anisimova M, Bielawski JP, Yang Z. Accuracy and power of bayes prediction of amino acid sites under positive selection. Mol Biol Evol. 2002;19(6):950–8.View ArticlePubMedGoogle Scholar
- Yang Z, Wong WS, Nielsen R. Bayes empirical bayes inference of amino acid sites under positive selection. Mol Biol Evol. 2005;22(4):1107–18.View ArticlePubMedGoogle Scholar
- Wilson DJ, Hernandez RD, Andolfatto P, Przeworski M. A population genetics-phylogenetics approach to inferring natural selection in coding sequences. PLoS Genet. 2011;7(12):e1002395.View ArticlePubMedPubMed CentralGoogle Scholar
- Wertheim JO, Murrell B, Smith MD, Kosakovsky Pond SL, Scheffler K. RELAX: detecting relaxed selection in a phylogenetic framework. Mol Biol Evol. 2015;32(3):820–32.View ArticlePubMedGoogle Scholar
- Wright S. Genetical structure of populations. Nature. 1950;166(4215):247–9.View ArticlePubMedGoogle Scholar
- Barreiro LB, Ben-Ali M, Quach H, Laval G, Patin E, Pickrell JK, et al. Evolutionary dynamics of human Toll-like receptors and their different contributions to host defense. PLoS Genet. 2009;5(7):e1000562.View ArticlePubMedPubMed CentralGoogle Scholar
- Voight BF, Kudaravalli S, Wen X, Pritchard JK. A map of recent positive selection in the human genome. PLoS Biol. 2006;4(3):e72.View ArticlePubMedPubMed CentralGoogle Scholar
- Watterson GA. On the number of segregating sites in genetical models without recombination. Theor Popul Biol. 1975;7(2):256–76.View ArticlePubMedGoogle Scholar
- Nei M, Li WH. Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc Natl Acad Sci U S A. 1979;76(10):5269–73.View ArticlePubMedPubMed CentralGoogle Scholar
- Mahajan A, Sim X, Ng HJ, Manning A, Rivas MA, Highland HM, et al. Identification and functional characterization of G6PC2 coding variants influencing glycemic traits define an effector transcript at the G6PC2-ABCB11 locus. PLoS Genet. 2015;11(1):e1004876.View ArticlePubMedPubMed CentralGoogle Scholar
- Wessel J, Chu AY, Willems SM, Wang S, Yaghootkar H, Brody JA, et al. Low-frequency and rare exome chip variants associate with fasting glucose and type 2 diabetes susceptibility. Nat Commun. 2015;6:5897.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang H, Liu L, Zhao J, Cui G, Chen C, Ding H, Wang DW. Large scale meta-analyses of fasting plasma glucose raising variants in GCK, GCKR, MTNR1B and G6PC2 and their impacts on type 2 diabetes mellitus risk. PLoS One. 2013;8(6):e67665.View ArticlePubMedPubMed CentralGoogle Scholar
- DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium, Asian Genetic Epidemiology Network Type 2 Diabetes (AGEN-T2D) Consortium, South Asian Type 2 Diabetes (SAT2D) Consortium, Mexican American Type 2 Diabetes (MAT2D) Consortium, Type 2 Diabetes Genetic Exploration by Nex-generation sequencing in muylti-Ethnic Samples (T2D-GENES) Consortium, Mahajan A, et al. Genome-wide trans-ancestry meta-analysis provides insight into the genetic architecture of type 2 diabetes susceptibility. Nat Genet. 2014;46(3):234–44.View ArticleGoogle Scholar
- Al-Nozha MM, Al-Maatouq MA, Al-Mazrou YY, Al-Harthi SS, Arafah MR, Khalil MZ, et al. Diabetes mellitus in Saudi Arabia. Saudi Med J. 2004;25(11):1603–10.PubMedGoogle Scholar
- Al-Othaimeen AI, Al-Nozha M, Osman AK. Obesity: an emerging problem in Saudi Arabia. Analysis of data from the National Nutrition Survey. East Mediterr Health J. 2007;13(2):441–8.PubMedGoogle Scholar
- Ionita-Laza I, Lee S, Makarov V, Buxbaum JD, Lin X. Sequence kernel association tests for the combined effect of rare and common variants. Am J Hum Genet. 2013;92(6):841–53.View ArticlePubMedPubMed CentralGoogle Scholar
- Echave J, Spielman SJ, Wilke CO. Causes of evolutionary rate variation among protein sites. Nat Rev Genet. 2016;17(2):109–21.View ArticlePubMedPubMed CentralGoogle Scholar
- Van de Peer Y, Maere S, Meyer A. The evolutionary significance of ancient genome duplications. Nat Rev Genet. 2009;10(10):725–32.View ArticlePubMedGoogle Scholar
- Smith JJ, Keinath MC. The sea lamprey meiotic map improves resolution of ancient vertebrate genome duplications. Genome Res. 2015;25(8):1081–90.View ArticlePubMedPubMed CentralGoogle Scholar
- Lucie M, Weiwei D, Stephane P, Sandrine SC. The five glucose-6-phosphatase paralogous genes are differentially regulated by insulin alone or combined with high level of amino acids and/or glucose in trout hepatocytes. Mol Biol Rep. 2016;43(4):207–11.View ArticlePubMedGoogle Scholar
- Prufer K, Racimo F, Patterson N, Jay F, Sankararaman S, Sawyer S, et al. The complete genome sequence of a Neanderthal from the Altai Mountains. Nature. 2014;505(7481):43–9.View ArticlePubMedGoogle Scholar
- Meyer M, Kircher M, Gansauge MT, Li H, Racimo F, Mallick S, et al. A high-coverage genome sequence from an archaic Denisovan individual. Science. 2012;338(6104):222–6.View ArticlePubMedPubMed CentralGoogle Scholar
- Prokopenko I, Langenberg C, Florez JC, Saxena R, Soranzo N, Thorleifsson G, et al. Variants in MTNR1B influence fasting glucose levels. Nat Genet. 2009;41(1):77–81.View ArticlePubMedGoogle Scholar
- Manning AK, Hivert MF, Scott RA, Grimsby JL, Bouatia-Naji N, Chen H, Rybin D, et al. A genome-wide approach accounting for body mass index identifies genetic variants influencing fasting glycemic traits and insulin resistance. Nat Genet. 2012;44(6):659–69.View ArticlePubMedPubMed CentralGoogle Scholar
- Hayes MG, Urbanek M, Hivert MF, Armstrong LL, Morrison J, Guo C, et al. Identification of HKDC1 and BACE2 as genes influencing glycemic traits during pregnancy through genome-wide association studies. Diabetes. 2013;62(9):3282–91.View ArticlePubMedPubMed CentralGoogle Scholar
- Hwang JY, Sim X, Wu Y, Liang J, Tabara Y, Hu C, et al. Genome-wide association meta-analysis identifies novel variants associated with fasting plasma glucose in East Asians. Diabetes. 2015;64(1):291–8.View ArticlePubMedGoogle Scholar
- Pare G, Chasman DI, Parker AN, Nathan DM, Miletich JP, Zee RY, Ridker PM. Novel association of HK1 with glycated hemoglobin in a non-diabetic population: a genome-wide evaluation of 14,618 participants in the Women’s Genome Health Study. PLoS Genet. 2008;4(12):e1000312.View ArticlePubMedPubMed CentralGoogle Scholar
- Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30(9):1312–3.View ArticlePubMedPubMed CentralGoogle Scholar
- Wernersson R, Pedersen AG. RevTrans: Multiple alignment of coding DNA from aligned amino acid sequences. Nucleic Acids Res. 2003;31(13):3537–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Pond SL, Frost SD, Muse SV. HyPhy: hypothesis testing using phylogenies. Bioinformatics. 2005;21(5):676–9.View ArticlePubMedGoogle Scholar
- Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59(3):307–21.View ArticlePubMedGoogle Scholar
- 1000 Genomes Project Consortium, Abecasis GR, Auton A, Brooks LD, DePristo MA, Durbin RM, Handsaker RE, Kang HM, Marth GT, McVean GA. An integrated map of genetic variation from 1,092 human genomes. Nature. 2012;491(7422):56–65.View ArticleGoogle Scholar
- Prado-Martinez J, Sudmant PH, Kidd JM, Li H, Kelley JL, Lorente-Galdos B, et al. Great ape genetic diversity and population history. Nature. 2013;499(7459):471–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Forni D, Cagliani R, Tresoldi C, Pozzoli U, De Gioia L, Filippi G, et al. An Evolutionary Analysis of Antigen Processing and Presentation across Different Timescales Reveals Pervasive Selection. PLoS Genet. 2014;10(3):e1004189.View ArticlePubMedPubMed CentralGoogle Scholar
- Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75.View ArticlePubMedPubMed CentralGoogle Scholar