BACKGROUND
the endoplasmic reticulum enzyme glucose-6-phosphatase catalyzes the hydrolysis of glucose-6-phosphate  to glucose and inorganic phosphate.
the enzyme is part of a multicomponent integral membrane system that includes the catalytic subunit  as well as transporters for glucose-6-phosphate, inorganic phosphate, and glucose  <cit> .
g6pase catalyzes the common terminal reaction in the gluconeogenic and glycogenolytic pathways, resulting in the release of glucose into the bloodstream  <cit> .
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 , usually referred to as the g6pc gene family  <cit> .
the protein products of the three genes display moderate sequence identity and a common topological organization with nine transmembrane domains and intralumenal catalytic residues  <cit> .
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  <cit> .
in humans, mutations in the gene cause glycogen storage disease type ia , 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  <cit> .
mutations in g6pc <dig> are also associated with pathology in humans.
thus, although the gene is ubiquitously expressed, its function is particularly important in white blood cells, and g6pc <dig> deficiency causes autosomal recessive severe congenital neutropenia type  <dig>   <cit> .
scn <dig> 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, g6pc <dig> only marginally contributes to the regulation of blood glucose levels or hepatic glycogen content  <cit> .
finally, g6pc <dig> is specifically expressed in pancreatic islets where its function is still incompletely understood  <cit> .
g6pc2
−/− mice display a reduction in blood glucose levels after a 6 h fast, whereas plasma insulin and glucagon concentrations are unaffected  <cit> .
these data led to the hypothesis that g6pc <dig> regulates the glycolytic flux by hydrolyzing g6p, thereby opposing the action of glucokinase.
g6pc <dig> and glucokinase are, therefore suggested to function as beta islet glucose sensors  <cit> .
in humans, common and rare variants in g6pc <dig> have been associated with fasting glucose  levels and with decreased insulin secretion during glucose tolerance tests .
this observation led to the suggestion that g6pc <dig> may also regulate the pulsatility of insulin secretion  <cit> .
variation in fg is clinically important in humans, as it is associated with the risk of developing type  <dig> diabetes  and ischemic heart disease  <cit>  as well as being an important determinant of offspring birth weight in pregnant women  <cit> .
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  <cit> .
mammals display a wide variety of diets, different lifestyles , and distinct energy requirements.
these characteristics influence the ability of a species to resist prolonged fasting  <cit> , 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  selection during mammalian evolution.
indeed, positive selection was previously demonstrated to act on genes with a role in carbohydrate absorption and digestion in mammals  <cit> .
in humans, aside from the textbook example of lactase persistence  <cit> , signals of diet-driven selection include variants in genes involved in starch and sucrose metabolism  <cit> , copy number variation at genes encoding salivary amylase   <cit> , as well as polymorphisms in genes that may be associated with the consumption of cooked food  <cit> .
in fact, humans likely underwent several dietary shifts associated with cultural innovations such as the use of fire for cooking   <cit> , the exploitation of starch-rich plant underground storage organs  <cit> , 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 g6pc <dig> variants in t2d susceptibility in a population with high incidence of metabolic disorders.
RESULTS
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  <dig> animal species through the ensembl compara utility  <cit>  indicated that a single g6pc gene is present in the drosophila genome, whereas lamprey  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.
overall, these observations suggest that the first duplication of an ancestral g6pc gene occurred during the vertebrate/invertebrate split and a second duplication took place either in the ancestor of all gnathostomata  or in the ancestor of bony vertebrates .
to more precisely map these duplication events, we constructed a phylogenetic tree using protein sequence information for the animal species included in the ensembl database plus additional organisms selected to resolve the timing of the duplication events .
results indicated that arthropods, mollusks, and echinoderms display one single g6pc gene, with the only exception of limulus polyphemus, which has two highly similar genes suggesting a recent duplication event in this lineage.
one g6pc gene is also observed in the hemichordate saccoglossus kowalevskii.
no g6pc gene was identified in the genomes of tunicates and cephalochordates, suggesting lineage-specific losses.fig.
<dig> maximum likelihood phylogenetic tree of metazoan g6pc proteins.
colored boxes indicate the class of each species , as reported in the legend.
black dots indicate bootstrap values greater than 50%
analysis of the g6pc phylogeny indicated that an initial duplication event in the lineage basal to all vertebrates originated g6pc <dig> and the g6pc/g6pc <dig> ancestor.
in lamprey, one of the two g6pc sequences clusters with g6pc <dig> proteins, whereas the other is basal to g6pc <dig> and g6pc , suggesting that the duplication events that originated g6pc and g6pc <dig> occurred after the split of gnathostomes and cyclostomes but before the divergence of cartilaginous and bony fishes, as the three callorhinchus milii sequences  indicate .
evolutionary analysis of the glucose-6-phosphatase  catalytic subunit gene family in mammals
we next analyzed in detail the evolutionary history of the three genes encoding g6pases in eutherian mammals.
to this aim, coding sequence information for ~ <dig> species were retrieved .
specifically, all available sequences with good coverage were retrieved for the study.
the rat sequence was not included for g6pc <dig>  as the gene is non functional in this rodent species  <cit> .
gard   <cit>  detected no recombination breakpoint in any alignment.
to obtain an estimate of the extent of functional constraint acting on these genes, we calculated the average non-synonymous substitution/synonymous substitution rate  using the single-likelihood ancestor counting  method  <cit> .
as is the case for most mammalian genes  <cit> , dn/ds was always lower than  <dig> , indicating that purifying selection is the major force shaping diversity at mammalian g6pase genes.
indeed, analysis based on the fixed effects likelihood  method  <cit>  detected a considerable proportion of negatively selected sites in all three genes .
the protein products of the three genes share a common topological structure, display considerable sequence identity, and perform the same molecular function, albeit in different cell types.
to test whether structural/functional constraints represent major drivers of molecular evolution, we used fel to calculate the normalized dn-ds value at each site and we correlated this parameter across corresponding sites .
although a significant correlation between dn-ds values was detected for g6pc and g6pc <dig> , as well as for g6pc and g6pc <dig> , the correlation coefficients were small .
no significant correlation was detected for the g6pc2-g6pc <dig> pair .table  <dig> average non-synonymous/synonymous substitution rate ratio  and percentage of negatively selected sites fot the three g6pase genes
g6pc
g6pt
g6pc2
igrp
g6pc3
ugrp
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,  <dig> independent gsd1a missense mutations involving  <dig> 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  than codons where no mutation has been described .
the same calculation was not performed for g6pc <dig> mutations, as too few of such mutations are actually known .
positive selection at the mammalian g6pc <dig> gene
positive selection may act on specific sites in a protein that is otherwise selectively constrained; to test for evidence of positive selection in the three g6pase genes, we applied likelihood ratio tests  implemented in the codeml program  <cit> .
the total tree length for eutherian mammals sequences varied between  <dig>  and  <dig>  ; these values are within an optimal accuracy range for codeml sites models  <cit> .
codeml was applied to compare models of gene evolution that allow  or disallow  a class of codons to evolve with dn/ds >  <dig>  as reported in table  <dig>  all null models were rejected in favor of the positive selection models for g6pc2; the same result was obtained using different codon frequency models  .
conversely, no evidence of positive selection was obtained for g6pc and g6pc <dig> .
these results indicate that g6pc <dig> alone evolved adaptively in mammals.
the bayes empirical bayes  analysis   <cit>  identified  <dig> codons showing strong evidence of positive selection ; most of these were also detected by fel or rel   <cit> .
with the exclusion of codon  <dig>  selected sites were located in the c-terminal portion of the protein, often within highly constrained regions .
human coding polymorphisms that modulate glycemic traits are mainly located in this c-terminal highly constrained region ; most of these variants affect codons that were targeted by negative selection during mammalian evolution .table  <dig> likelihood ratio test statistics for models of variable selective pressure among sites in g6pc2
am1a is a nearly neutral model that assumes one ω class between  <dig> and  <dig> and one class with ω = 1; m2a  is the same as m1a plus an extra class of ω > 1
bm <dig> is a null model that assumes that 0 < ω <  <dig> is beta distributed among sites; m <dig>  is the same as m <dig> but also includes an extra category of sites with ω > 1
cm8a is the same as m <dig>  except that the 11th category cannot allow positive selection, but only neutral evolution
d2Δlnl: twice the difference of the natural logs of the maximum likelihood of the models being compared
fig.
<dig> evolutionary analysis of g6pase genes.
a g6pc <dig> is shown with its predicted membrane topology; protein regions are coloured in hues of blue according to the percentage of negatively selected sites .
positively selected sites in the mammalian phylogeny  and in homininae  are reported on the structure.
missense variants associated with fg are shown in red.
asterisks denote negatively selected sites.
the glycosylation site is also shown.
b violin plots of selection coefficients  for the three g6pase genes.
selection coefficients  are classified as strongly beneficial , moderately beneficial , weakly beneficial , neutral , weakly deleterious , moderately deleterious , strongly deleterious , and inviable .
c phylogenetic tree for primate g6pc <dig> genes.
branches are color-coded according to relax results: blue, significant evidence of relaxed selection; orange, no significant evidence of relaxation
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  <cit>  that jointly uses intra-specific variation and inter-specific diversity to estimate the distribution of fitness effects  along coding regions.
gammamap categorizes codons into  <dig> classes of γ, ranging from strongly beneficial  to inviable ; a γ equal to  <dig> indicates neutrality.
the overall distribution of selection coefficients indicated that g6pc evolved under strong purifying selection in all lineages .
this was also the case for g6pc <dig> in non-human primates , whereas the human gene showed weaker constraint .
finally, the distribution of fitness effects for g6pc <dig> 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  <dig>  in contrast, the gorilla lineage showed evidence of strong purifying selection .
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  <cit>  to the g6pc <dig> primate phylogeny .
results were consistent with relaxed selection on the human/chimpanzee branches , but not on the gorilla lineage  .
the same analysis for the human g6pc <dig> branch revealed no relaxation of selective pressure .
gammamap also identified two positively selected codons  for human g6pc <dig> .
one selected codon was also identified for human g6pc <dig> , 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  <dig> genomes phase  <dig> data for yoruba, european, and chinese we calculated pairwise fst  <cit> , an estimate of population genetic differentiation.
we also performed the dind  and ihs  tests  <cit>  for all snps mapping to these genes.
statistical significance  for the fst statistic and for the dind test was obtained by deriving empirical distributions.
for the ihs test, absolute values higher than  <dig> were considered as significant  <cit> .
no snp in any g6pase gene reached statistical significance  for both fst and for the dind tests, and none had an |ihs| higher than  <dig>  overall, these results indicate that no variant/haplotype can be confidently called as positively selected.
likewise, nucleotide diversity  for the entire gene regions was unexceptional if compared to those calculated for a reference set of  <dig> genes.
we conclude that g6pase genes did not represent selection targets in recent human history.
association of g6pc <dig> variants with t2d
several genome-wide association studies  have identified a functional non-coding variant  in g6pc <dig> that is associated with fasting glucose  levels .
more recently, multiple rare and common coding variants in this gene were shown to influence fg  <cit> .
as mentioned above, all these coding variants are located in the two terminal exons of g6pc <dig>  where most sites that are positively selected in mammals also map .
the best characterized variants  exert an effect independent of each other and of the gwas snp, indicating that haplotype analysis rather than single variant association is better suited to assess the contribution of g6pc <dig> variants to metabolic traits  <cit> .
despite their replicated effect on fg, the contribution of rare and common g6pc <dig> variants to t2d susceptibility has remained controversial  <cit> .
we thus investigated a possible role for g6pc <dig> variants in modulating the susceptibility to t2d in subjects from saudi arabia, a region with a high prevalence of metabolic disorders, including t2d  <cit> .
specifically, we resequenced the two terminal exons of g6pc <dig>  in  <dig> subjects from saudi arabia,  <dig> of whom suffering from t2d .
to limit phenotype heterogeneity only non-obese individuals  were included.
the rs <dig> gwas variant was also genotyped.fig.
<dig> linkage disequilibrium  plots.
the ld plot was constructed with haploview  <dig>  and displays r
<dig> values  for the polymorphic variants we identified.
ld plots of the common variants for ceu, yri and chb is also shown.
the exon-intron structure of g6pc <dig>  is also shown together with the two regions we resequenced
no novel missense or nonsense variant was detected in either t2d subjects or healthy controls  and the frequency of known rare missense and nonsense variants was not significantly different in t2d and hc .
two common missense variants were nevertheless detected in the last g6pc <dig> exon: rs <dig>  and rs <dig> .
the two variants display very limited linkage disequilibrium  .
to address their contribution to t2d risk, logistic regression using age, sex, and bmi as covariates were used.
after fdr correction for multiple tests, no association with t2d was observed ; conditioning on the gwas variant, though, revealed a significant association for the two missense variants .
haplotype analysis using the same covariates indicated above detected two haplotypes significantly associated with t2d .
both the predisposing and the protective haplotype carry the glucose-raising allele at rs <dig>  the predisposing haplotype also includes the loss-of-function l <dig> allele  and the minor allele  at rs <dig> .
these results should be regarded as preliminary due to the small sample size.table  <dig> association of g6pc <dig> variants with t2d
g6pc <dig> haplotype analysis
finally, to assess the effect of rare and common g6pc <dig> variants on t2d risk, we applied a snp-set based method, the sequence kernel association test   <cit> .
skat was run either by inclusion of all variants identified through re-sequencing  or by limiting analysis to missense snps plus the gwas variant .
no significant association was detected.
however, as for single-variant associations, the power of skat is limited when small samples are analyzed  <cit> .
discussion
in this study we have analyzed the evolutionary history of three genes  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  <cit> .
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  <cit> .
in line with this view, we found all g6pase genes to display an overall dn/ds lower than  <dig>  indicating a preponderance of negative selection.
recent evidence showed that structural and folding requirements  represent major determinants of the evolutionary rate at protein sites  <cit> .
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  <cit> .
nonetheless, we reasoned that because the three proteins share considerable identity in terms of amino acid sequence and the same topological organization  <cit> , 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  <cit>  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  <cit> .
it is thus possible that g6pc <dig> and the g6pc/g6pc <dig> ancestor originated and were fixed after whole genome duplication in the ancestral vertebrate.
however, the basal position of one lamprey sequence with respect to gnathostome g6pc and g6pc <dig> 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  <dig> g6pc paralogs, possibly as a results of a whole genome duplication that occurred in the ancestor of teleosts  <cit> .
a similar observation was reported for the rainbow trout, a glucose-intolerant fish, which displays  <dig> g6pc genes possibly fixed in this species after the salmonid-specific whole genome duplication  <cit> .
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 g6pc <dig> 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 g6pc <dig> 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 g6pc <dig> 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 c <dig> and a <dig> variant were already present in the genomes of neandertals and denisovans  <cit> .
these observations suggest that, as for other variants in metabolic genes  <cit> , 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 g6pc <dig> 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  <dig> nonsynonymous substitutions , including the positively selected  <dig> site, in the absence of synonymous substitution.
three of these changes are clustered in ~ <dig> amino acid region  suggesting that, for unknown reasons, this protein portion is tolerant to change in humans.
to date, no snc <dig> missense mutation has been described at these sites.
among the three g6pase genes, mutations in g6pc <dig> 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 g6pc <dig> is a pseudogene in rats.
moreover, recent functional data indicated that coding variants that reduce the expression of g6pc <dig>  most likely by impairing its folding, segregate at appreciable frequency in human populations  <cit> .
notably, variants in g6pc <dig> have been consistently associated with fg levels, whereas their contribution to t2d risk remains controversial.
in particular, the rs <dig> snp is one of the strongest signals associated to fg , and one of the most commonly replicated in large-scale analyses .
moreover, the variant was shown to be functional and to modulate g6pc <dig> pre-mrna splicing  <cit> .
although this latter finding does not necessarily imply that rs <dig> is the causal variant, the effect of the glucose-raising allele  on increased splicing efficiency is suggestive  <cit> .
however, distinct studies found either no association of rs <dig> with t2d risk  <cit>  or indicated a weak protective effect of the glucose-increasing allele  <cit> .
recently, mahajan and coworkers reported a glucose-increasing effect of the common v <dig>  allele that modestly increases the risk of t2d as well  <cit> .
the authors suggested that association analysis for g6pc <dig> should be performed through haplotype reconstruction as multiple rare and common variants independently affect fg levels, and the direction of effect for rs <dig> is reversed when analysis is conditioned on rs <dig>  <cit> .
nonetheless, most large-scale analyses of t2d susceptibility performed single variant association tests, rather than haplotype inference, leaving the role of g6pc <dig> 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 rs <dig>  variant, that is not covered in exome chip arrays and was thus not analyzed in recent association studies of g6pc <dig> variants for fg levels  <cit> .
in a genome-wide meta-analysis  <cit> , rs <dig> showed a strong association with fg , which is likely independent of the lead variant rs <dig>  as their ld is low in all populations  .
the functional effect of the s342c substitution is presently unknown.
codon  <dig> 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 .
these observations suggest that the derived s <dig> allele impairs g6pc <dig> function.
surprisingly, though, the v <dig> allele which also involves a negatively selected site and represents the ancestral state conserved in all mammals , was recently shown to result in reduced function  <cit> .
indeed, g6pc <dig> molecules carrying the v <dig> allele are expressed at lower abundance due to proteasomal degradation  <cit> .
this observation indicates that the functional effect of g6pc <dig> variants is difficult to predict, and in the case of the s <dig> substitution will need experimental testing.
the data we report herein, although preliminary, may help reconcile the contrasting results obtained for rs <dig> on t2d risk, as its effect might depend on haplotype context and may vary in different populations depending on ld between rs <dig> and other functional variants.
clearly, further studies will be necessary to confirm the role of g6pc <dig> 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 g6pc <dig>  and in the intergenic spacer downstream the transcription end site of the gene  were also associated with fg  <cit> .
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 g6pc <dig> genetic variability to t2d risk.
CONCLUSIONS
in conclusion, we detected pervasive positive selection at mammalian g6pc <dig> genes, with almost all selected sites located in the c-terminal portion of the protein.
we then investigated a possible role for g6pc <dig> 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 g6pc <dig> haplotypes modulate susceptibility to t2d.

 the endoplasmic reticulum enzyme glucose-6-phosphatase catalyzes the hydrolysis of glucose-6-phosphate to glucose and inorganic phosphate. the enzyme is part of a multicomponent integral membrane system that includes the catalytic subunit as well as transporters for glucose-6-phosphate, inorganic phosphate, and glucose. 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, usually referred to as the g6pc gene family. 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. in humans, mutations in the gene cause glycogen storage disease type ia, 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. mutations in g6pc are also associated with pathology in humans. thus, although the gene is ubiquitously expressed, its function is particularly important in white blood cells, and g6pc deficiency causes autosomal recessive severe congenital neutropenia type. scn 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, g6pc only marginally contributes to the regulation of blood glucose levels or hepatic glycogen content. finally, g6pc is specifically expressed in pancreatic islets where its function is still incompletely understood. g6pc2 −/− mice display a reduction in blood glucose levels after a 6 h fast, whereas plasma insulin and glucagon concentrations are unaffected. these data led to the hypothesis that g6pc regulates the glycolytic flux by hydrolyzing g6p, thereby opposing the action of glucokinase. g6pc and glucokinase are, therefore suggested to function as beta islet glucose sensors. in humans, common and rare variants in g6pc have been associated with fasting glucose levels and with decreased insulin secretion during glucose tolerance tests. this observation led to the suggestion that g6pc may also regulate the pulsatility of insulin secretion. variation in fg is clinically important in humans, as it is associated with the risk of developing type diabetes and ischemic heart disease 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, 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 selection during mammalian evolution. indeed, positive selection was previously demonstrated to act on genes with a role in carbohydrate absorption and digestion in mammals. 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, copy number variation at genes encoding salivary amylase, 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, 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 g6pc 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 animal species through the ensembl compara utility indicated that a single g6pc gene is present in the drosophila genome, whereas lamprey 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. overall, these observations suggest that the first duplication of an ancestral g6pc gene occurred during the vertebrate/invertebrate split and a second duplication took place either in the ancestor of all gnathostomata or in the ancestor of bony vertebrates. to more precisely map these duplication events, we constructed a phylogenetic tree using protein sequence information for the animal species included in the ensembl database plus additional organisms selected to resolve the timing of the duplication events. results indicated that arthropods, mollusks, and echinoderms display one single g6pc gene, with the only exception of limulus polyphemus, which has two highly similar genes suggesting a recent duplication event in this lineage. one g6pc gene is also observed in the hemichordate saccoglossus kowalevskii. no g6pc gene was identified in the genomes of tunicates and cephalochordates, suggesting lineage-specific losses.fig. maximum likelihood phylogenetic tree of metazoan g6pc proteins. colored boxes indicate the class of each species, as reported in the legend. black dots indicate bootstrap values greater than 50% analysis of the g6pc phylogeny indicated that an initial duplication event in the lineage basal to all vertebrates originated g6pc and the g6pc/g6pc ancestor. in lamprey, one of the two g6pc sequences clusters with g6pc proteins, whereas the other is basal to g6pc and g6pc, suggesting that the duplication events that originated g6pc and g6pc occurred after the split of gnathostomes and cyclostomes but before the divergence of cartilaginous and bony fishes, as the three callorhinchus milii sequences indicate. evolutionary analysis of the glucose-6-phosphatase catalytic subunit gene family in mammals we next analyzed in detail the evolutionary history of the three genes encoding g6pases in eutherian mammals. to this aim, coding sequence information for ~ species were retrieved. specifically, all available sequences with good coverage were retrieved for the study. the rat sequence was not included for g6pc as the gene is non functional in this rodent species. gard detected no recombination breakpoint in any alignment. to obtain an estimate of the extent of functional constraint acting on these genes, we calculated the average non-synonymous substitution/synonymous substitution rate using the single-likelihood ancestor counting method. as is the case for most mammalian genes, dn/ds was always lower than, indicating that purifying selection is the major force shaping diversity at mammalian g6pase genes. indeed, analysis based on the fixed effects likelihood method detected a considerable proportion of negatively selected sites in all three genes. the protein products of the three genes share a common topological structure, display considerable sequence identity, and perform the same molecular function, albeit in different cell types. to test whether structural/functional constraints represent major drivers of molecular evolution, we used fel to calculate the normalized dn-ds value at each site and we correlated this parameter across corresponding sites. although a significant correlation between dn-ds values was detected for g6pc and g6pc, as well as for g6pc and g6pc, the correlation coefficients were small. no significant correlation was detected for the g6pc2-g6pc pair.table average non-synonymous/synonymous substitution rate ratio and percentage of negatively selected sites fot the three g6pase genes g6pc g6pt g6pc2 igrp g6pc3 ugrp 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, independent gsd1a missense mutations involving 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 than codons where no mutation has been described. the same calculation was not performed for g6pc mutations, as too few of such mutations are actually known. positive selection at the mammalian g6pc gene positive selection may act on specific sites in a protein that is otherwise selectively constrained; to test for evidence of positive selection in the three g6pase genes, we applied likelihood ratio tests implemented in the codeml program. the total tree length for eutherian mammals sequences varied between and; these values are within an optimal accuracy range for codeml sites models. codeml was applied to compare models of gene evolution that allow or disallow a class of codons to evolve with dn/ds >  as reported in table  all null models were rejected in favor of the positive selection models for g6pc2; the same result was obtained using different codon frequency models. conversely, no evidence of positive selection was obtained for g6pc and g6pc. these results indicate that g6pc alone evolved adaptively in mammals. the bayes empirical bayes analysis identified codons showing strong evidence of positive selection; most of these were also detected by fel or rel. with the exclusion of codon selected sites were located in the c-terminal portion of the protein, often within highly constrained regions. human coding polymorphisms that modulate glycemic traits are mainly located in this c-terminal highly constrained region; most of these variants affect codons that were targeted by negative selection during mammalian evolution.table likelihood ratio test statistics for models of variable selective pressure among sites in g6pc2 am1a is a nearly neutral model that assumes one ω class between and and one class with ω = 1; m2a is the same as m1a plus an extra class of ω > 1 bm is a null model that assumes that 0 < ω <  is beta distributed among sites; m is the same as m but also includes an extra category of sites with ω > 1 cm8a is the same as m except that the 11th category cannot allow positive selection, but only neutral evolution d2Δlnl: twice the difference of the natural logs of the maximum likelihood of the models being compared fig. evolutionary analysis of g6pase genes. a g6pc is shown with its predicted membrane topology; protein regions are coloured in hues of blue according to the percentage of negatively selected sites. positively selected sites in the mammalian phylogeny and in homininae are reported on the structure. missense variants associated with fg are shown in red. asterisks denote negatively selected sites. the glycosylation site is also shown. b violin plots of selection coefficients for the three g6pase genes. selection coefficients are classified as strongly beneficial, moderately beneficial, weakly beneficial, neutral, weakly deleterious, moderately deleterious, strongly deleterious, and inviable. c phylogenetic tree for primate g6pc genes. branches are color-coded according to relax results: blue, significant evidence of relaxed selection; orange, no significant evidence of relaxation 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 along coding regions. gammamap categorizes codons into classes of γ, ranging from strongly beneficial to inviable; a γ equal to indicates neutrality. the overall distribution of selection coefficients indicated that g6pc evolved under strong purifying selection in all lineages. this was also the case for g6pc in non-human primates, whereas the human gene showed weaker constraint. finally, the distribution of fitness effects for g6pc 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 in contrast, the gorilla lineage showed evidence of strong purifying selection. 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 g6pc primate phylogeny. results were consistent with relaxed selection on the human/chimpanzee branches, but not on the gorilla lineage. the same analysis for the human g6pc branch revealed no relaxation of selective pressure. gammamap also identified two positively selected codons for human g6pc. one selected codon was also identified for human g6pc, 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 genomes phase data for yoruba, european, and chinese we calculated pairwise fst, an estimate of population genetic differentiation. we also performed the dind and ihs tests for all snps mapping to these genes. statistical significance for the fst statistic and for the dind test was obtained by deriving empirical distributions. for the ihs test, absolute values higher than were considered as significant. no snp in any g6pase gene reached statistical significance for both fst and for the dind tests, and none had an |ihs| higher than overall, these results indicate that no variant/haplotype can be confidently called as positively selected. likewise, nucleotide diversity for the entire gene regions was unexceptional if compared to those calculated for a reference set of genes. we conclude that g6pase genes did not represent selection targets in recent human history. association of g6pc variants with t2d several genome-wide association studies have identified a functional non-coding variant in g6pc that is associated with fasting glucose levels. more recently, multiple rare and common coding variants in this gene were shown to influence fg. as mentioned above, all these coding variants are located in the two terminal exons of g6pc where most sites that are positively selected in mammals also map. the best characterized variants exert an effect independent of each other and of the gwas snp, indicating that haplotype analysis rather than single variant association is better suited to assess the contribution of g6pc variants to metabolic traits. despite their replicated effect on fg, the contribution of rare and common g6pc variants to t2d susceptibility has remained controversial. we thus investigated a possible role for g6pc variants in modulating the susceptibility to t2d in subjects from saudi arabia, a region with a high prevalence of metabolic disorders, including t2d. specifically, we resequenced the two terminal exons of g6pc in subjects from saudi arabia, of whom suffering from t2d. to limit phenotype heterogeneity only non-obese individuals were included. the rs gwas variant was also genotyped.fig. linkage disequilibrium plots. the ld plot was constructed with haploview and displays r values for the polymorphic variants we identified. ld plots of the common variants for ceu, yri and chb is also shown. the exon-intron structure of g6pc is also shown together with the two regions we resequenced no novel missense or nonsense variant was detected in either t2d subjects or healthy controls and the frequency of known rare missense and nonsense variants was not significantly different in t2d and hc. two common missense variants were nevertheless detected in the last g6pc exon: rs and rs. the two variants display very limited linkage disequilibrium. to address their contribution to t2d risk, logistic regression using age, sex, and bmi as covariates were used. after fdr correction for multiple tests, no association with t2d was observed; conditioning on the gwas variant, though, revealed a significant association for the two missense variants. haplotype analysis using the same covariates indicated above detected two haplotypes significantly associated with t2d. both the predisposing and the protective haplotype carry the glucose-raising allele at rs the predisposing haplotype also includes the loss-of-function l allele and the minor allele at rs. these results should be regarded as preliminary due to the small sample size.table association of g6pc variants with t2d g6pc haplotype analysis finally, to assess the effect of rare and common g6pc variants on t2d risk, we applied a snp-set based method, the sequence kernel association test. skat was run either by inclusion of all variants identified through re-sequencing or by limiting analysis to missense snps plus the gwas variant. no significant association was detected. however, as for single-variant associations, the power of skat is limited when small samples are analyzed. discussion in this study we have analyzed the evolutionary history of three genes 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 indicating a preponderance of negative selection. recent evidence showed that structural and folding requirements 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 g6pc and the g6pc/g6pc ancestor originated and were fixed after whole genome duplication in the ancestral vertebrate. however, the basal position of one lamprey sequence with respect to gnathostome g6pc and g6pc 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 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 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 g6pc 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 g6pc 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 g6pc 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 c and a variant were already present in the genomes of neandertals and denisovans. 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 g6pc 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 nonsynonymous substitutions, including the positively selected site, in the absence of synonymous substitution. three of these changes are clustered in ~ amino acid region suggesting that, for unknown reasons, this protein portion is tolerant to change in humans. to date, no snc missense mutation has been described at these sites. among the three g6pase genes, mutations in g6pc 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 g6pc is a pseudogene in rats. moreover, recent functional data indicated that coding variants that reduce the expression of g6pc most likely by impairing its folding, segregate at appreciable frequency in human populations. notably, variants in g6pc have been consistently associated with fg levels, whereas their contribution to t2d risk remains controversial. in particular, the rs snp is one of the strongest signals associated to fg, and one of the most commonly replicated in large-scale analyses. moreover, the variant was shown to be functional and to modulate g6pc pre-mrna splicing. although this latter finding does not necessarily imply that rs is the causal variant, the effect of the glucose-raising allele on increased splicing efficiency is suggestive. however, distinct studies found either no association of rs with t2d risk or indicated a weak protective effect of the glucose-increasing allele. recently, mahajan and coworkers reported a glucose-increasing effect of the common v allele that modestly increases the risk of t2d as well. the authors suggested that association analysis for g6pc should be performed through haplotype reconstruction as multiple rare and common variants independently affect fg levels, and the direction of effect for rs is reversed when analysis is conditioned on rs. nonetheless, most large-scale analyses of t2d susceptibility performed single variant association tests, rather than haplotype inference, leaving the role of g6pc 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 rs variant, that is not covered in exome chip arrays and was thus not analyzed in recent association studies of g6pc variants for fg levels. in a genome-wide meta-analysis, rs showed a strong association with fg, which is likely independent of the lead variant rs as their ld is low in all populations. the functional effect of the s342c substitution is presently unknown. codon 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. these observations suggest that the derived s allele impairs g6pc function. surprisingly, though, the v allele which also involves a negatively selected site and represents the ancestral state conserved in all mammals, was recently shown to result in reduced function. indeed, g6pc molecules carrying the v allele are expressed at lower abundance due to proteasomal degradation. this observation indicates that the functional effect of g6pc variants is difficult to predict, and in the case of the s substitution will need experimental testing. the data we report herein, although preliminary, may help reconcile the contrasting results obtained for rs on t2d risk, as its effect might depend on haplotype context and may vary in different populations depending on ld between rs and other functional variants. clearly, further studies will be necessary to confirm the role of g6pc 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 g6pc and in the intergenic spacer downstream the transcription end site of the gene were also associated with fg. 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 g6pc genetic variability to t2d risk.  in conclusion, we detected pervasive positive selection at mammalian g6pc genes, with almost all selected sites located in the c-terminal portion of the protein. we then investigated a possible role for g6pc 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 g6pc haplotypes modulate susceptibility to t2d. 
nonetheless, site-wise estimation of evolutionary rates at corresponding sites revealed weak correlations, suggesting that mammalian g6pases have evolved different structural features over time. association analysis with two common missense variants revealed a weak but significant association for both snps when analyses were conditioned on rs previously identified in a gwas for fg. 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. we also detected pervasive positive selection at mammalian g6pc 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 ~ subjects from saudi arabia, of whom suffering from type diabetes. we analyzed the evolutionary history of g6pase genes. in most mammals, different g6pc subunits are encoded by three paralogous genes. the frequency of rare missense and nonsense variants was not significantly different in t2d and controls. 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. we detected pervasive positive selection at mammalian g6pc genes and we suggest that distinct haplotypes at the g6pc locus modulate susceptibility to t2d. g6pc2g6pcg6pc3natural selectionassociation analysistype diabetesissue-copyright-statement© the author 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. mutations in g6pc and g6pc are responsible for human mendelian diseases, whereas variants in g6pc are associated with fasting glucose 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 g6pc 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 ~ subjects from saudi arabia, of whom suffering from type diabetes. the frequency of rare missense and nonsense variants was not significantly different in t2d and controls. association analysis with two common missense variants revealed a weak but significant association for both snps when analyses were conditioned on rs 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 g6pc genes and we suggest that distinct haplotypes at the g6pc locus modulate susceptibility to t2d. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords g6pc2g6pcg6pc3natural selectionassociation analysistype diabetesissue-copyright-statement© the author 2017 