Comparison of the amino acid sequences and syntenic analyses of sarcopterygian g6pc1

The CLUSTAL W program embedded in MEGA (version 7.0.21) (Kumar et al., 2016) was used with default settings to align the nucleotide and deduced amino acid sequences of the sarcopterygian g6pc1. Syntenic analyses of sarcopterygian g6pc1 were performed using the genome data deposited in GenBank and Ensembl (see  Supplementary Table S1 (zs210113_TableS1.xlsx)).

Construction of molecular phylogenetic tree of vertebrate g6pc1

Molecular phylogenetic trees of vertebrate g6pc1 were constructed together with the g6pc2 and g6pc3 genes. In all trees, g6pc3 genes were used to root the trees. This is because the vertebrate g6pc family first bifurcates between g6pc1/g6pc2 and g6pc3 (Al-Daghri et al., 2017). Codon alignments of the CDS nucleotide sequences were constructed using the CLUSTAL W program embedded in MEGA (version 7.0.21) (Kumar et al., 2016) and were subsequently used for tree construction.

Three trees were selected for the study. First, a neighbor-joining (NJ) tree was drawn based on amino acid sequences to determine whether the results of Marandel et al. (2017) could be reproduced. A nucleotide-based tree was then constructed using the maximum likelihood (ML) method to re-examine the evolutionary history of vertebrate g6pc1. These trees contain both actinopterygian and sarcopterygian genes. To infer the evolution of sarcopterygian g6pc1 genes, another ML tree was drawn based exclusively on the nucleotide sequences of sarcopterygian genes. The amino acid-based NJ tree was drawn by selecting the “partial deletion” option for gap treatment with 95% of site coverage cutoff. Nucleotide-based trees were drawn by selecting the “partial deletion” option for gap treatment with 95% of site coverage cutoff. A total of 340 positions were included in the actinopterygian-sarcopterygian trees, whereas 342 positions were included in the sarcopterygian ML tree. Before constructing ML trees, “Find Best DNA/Protein Models” was run for the selection of “Substitution Model” and “Rates and Patterns”. The best setting for the actinopterygian + sarcopterygian nucleotide-based tree was the General Time Reversible model for “Substitution Model” and Gamma Distributed with Invariant sites (G + I) for “Rates and Patterns” (GTR + G + I). The best setting for the sarcopterygian nucleotide-based tree was the General Time Reversible model for the “Substitution Model” and Gamma Distributed (G) for “Rates and Patterns” (GTR + G). The number of bootstrap replicates for the estimation of the robustness of the internal nodes was set to 1000 for NJ and 100 for ML trees. The full species names and GenBank accession numbers are summarized in  Supplementary Table S1 (zs210113_TableS1.xlsx).

Comparison of the nucleotide sequences upstream of g6pc1

Human and mouse g6pc1 genes have multiple response elements in the region within –1000 base pairs (bp) from the transcription start site (Schmoll et al., 1996; Streeper et al., 1997, 1998; Qiu et al., 2017). We predicted that the genomic environment corresponding to this region would have a significant effect on transcriptional regulation. In this study, nucleotide sequences –1000 to –1 bp upstream of the transcription start site were aligned using the LAGAN program (Brudno et al., 2003) in the Vista package ( http://genome.lbl.gov/vista/mvista/submit.shtml) (Frazer et al., 2004).

Animals

Fasting experiments

Fourteen male Japanese geckos were divided into two groups of seven individuals each. After the acclimation period, one group was fasted for 23 days, whereas the other group was fed ad libitum.

Comparison of seasonal differences in free-living population

Male and female Japanese geckos were captured on 25 April 2017 (n = 5 for each sex) and on 26 September 2017 (n = 6 for each sex). All of the geckos were sacrificed on the same day that they were captured.

Administration of 17β-estradiol (E2)

Gecko species can store sperm in the oviduct for months after copulation (Murphy-Walker and Haley, 1996; Yamamoto and Ota, 2006). Therefore, the inadvertent onset of vitellogenesis and concomitant increase in blood estrogen levels are possible. In this experiment, such an effect was avoided by using males. Sixteen male Japanese geckos were divided into two groups of eight individuals each. After the acclimation period, the E2 treatment group was intraperitoneally injected with 17β-estradiol in sesame oil at a dose of 1 µg/g body weight every alternate day for 2 weeks. The control group was intraperitoneally injected with sesame oil alone. Both groups were provided with food and water ad libitum. Two individuals in the E2 treatment group died at 2 and 13 days after the first administration, respectively. Therefore, subsequent analyses were performed on six specimens from the E2 treatment group and eight specimens from the control group.

RNA extraction and cDNA synthesis

Total RNA was extracted from the liver samples using ISOGEN (Nippon Gene, Tokyo, Japan). For the construction of the DNA complementary to RNA (cDNA) templates, 3 µg of total RNA was used. The cDNA templates for quantitative real-time PCR were constructed as follows: denatured total RNA was reverse-transcribed using an oligo (oligodeoxyribonucleotide) (dT) primer and M-MLV reverse transcriptase (Promega, Madison, WI). The reaction volume contained oligo (dT) primers at 5 pM, 2 mM of each deoxyribonucleoside triphosphate (dNTP), 200 units of M-MLV reverse transcriptase, and 1 × M-MLV buffer (Promega). The volume was brought to 20 µL with nuclease-free water (Qiagen, Tokyo, Japan). The reaction was performed at 42°C for 90 min, and then 70°C for 15 min.

Quantitative real-time PCR

Five microliters of each 25- or 100-fold-diluted cDNA from liver samples was amplified in a 20 µl reaction mixture containing each primer (0.3 µM) and the KAPA SYBR Fast qPCR Kit (KAPA Biosystems, Century City, South Africa) using a LightCycler 480 Instrument II (Roche Diagnostic). The primers were designed to span at least one intron of more than 400 bp based on the genome sequences of the Japanese gecko in the NCBI database. The primer sequences used are listed in Table 1. The PCR reaction conditions were as follows: 95°C for 5 min, 45 cycles of 95°C for 10 s, 56°C for 10 s, and 72°C for 10 s, followed by melting curve analysis by which the PCR product was confirmed. The PCR product was also confirmed using RNA samples without RT reaction as a template for PCR. The specificity of the primers for each g6pc1 gene was confirmed by sequencing PCR products. The reaction efficiency was greater than 90% for all the primer sets. The amplicons from each gene were designed as follows: 115 bp for g6pc1-1; 142 bp for g6pc1-2; 112 bp for pck2; 88 bp for 36b4; 121 bp for rpl8.

Table 1.

Primers used for quantitative RT-PCR analyses.

Fig. 1.

Amino acid sequences deduced from sarcopterygian g6pc1 genes. The amino acid sequences of G6PC1, deduced from sarcopterygian genes, were aligned. It should be noted that only functionally important regions are presented. Species were selected from the major extant sarcopterygian lineages (Eutherian mammal, Marsupial, Squamata, Aves, Tuatara, Crocodilia, Chelonia, Amphibia, and coelacanth). Amino acids corresponding to active centers for enzymatic activities (Lys76, Arg83, His119, Arg170, and His176 in human G6PC1) are denoted by open triangles. All species and paralogs conserved both the phosphatase motif and the active centers.

Reference genes for normalization were selected using the NormFinder (Version 0.953) algorithm (Andersen et al., 2004), which is available as an Excel add-in at  https://moma.dk/normfindersoftware. 36b4, ef1a1, and actb levels were examined in the fasting experiment. Among these genes, 36b4 was selected as the most stable gene, with a stable value of 0.013. To compare the free-living individuals and for the E2 administration experiment, 36b4, ef1a1, actb, sdha, and rpl8 were examined. In both experiments, the geometric means of 36b4 and rpl8 were selected as the most stable reference, with stable values of 0.006 (free-living population) and 0.001 (E2 administration).

Search for hormone responsive elements upstream of Japanese gecko g6pc1s

Two types of hormone responsive elements, estrogen responsive element (ERE) and progesterone responsive element (PRE), were searched in the genomic regions upstream of the Japanese gecko g6pc1 genes. Genomic sequence information corresponding to 1 kilobase pair (kbp) upstream of the transcription start site was obtained from the GenBank database, and ERE and PRE half sites (AGGTCA and TGACCT for the former, AGAACA and TGTTCT for the latter) were manually searched with the aid of Windows Notepad.

Statistical analysis

All data are expressed as the mean ± standard error of the mean (SEM). The differences of the expression levels of the genes between the groups were analyzed using the Mann–Whitney U test. All analyses were conducted using KyPlot 5.0 software (KyensLab, Tokyo, Japan).

Comparison of deduced amino acid sequences and syntenic analyses of g6pc1 genes

The amino acid sequences of g6pc1 deduced from the genome data revealed that all sarcopterygian paralogs possessed a phosphatase motif. They also conserved five active center sites corresponding to Lys76, Arg83, His119, Arg170, and His176 of human G6PC1 (Chou and Mansfield, 2008) (Fig. 1).

Syntenic analysis of g6pc1 revealed that the coelacanth had a single g6pc1 gene (Fig. 2). According to GenBank annotation, the western clawed frog had two g6pc1 and one g6pc2 genes arranged in tandem on chromosome 10. However, this “g6pc2” gene was not included in the g6pc2 clade in the phylogenetic trees constructed in this study (Figs. 3, 4). In addition, no other sarcopterygian species had the g6pc1 and g6pc2 genes on the same chromosome (Fig. 2). Therefore, we designated this “g6pc2” as a g6pc1 gene. Among amniotes, mammals and aves had one g6pc1 gene, whereas reptiles had two or three genes. Among reptiles, squamates, crocodilians, and the softshell turtle had two genes. In contrast, the tuatara, green sea turtle, and western painted turtle had three. In all species investigated, the regions containing g6pc1 genes were flanked by the aoc and aarsd1 genes (Fig. 2).

Molecular phylogenetic analyses

To infer the phylogenetic relationships among vertebrate g6pc1 paralogs, three molecular phylogenetic trees were constructed. First, an NJ tree was constructed based on the amino acid sequences of the vertebrate g6pc family genes. This tree includes both actinopterygian and sarcopterygian species. In this NJ tree, the gnathostome g6pc1 paralogs formed a monophyletic clade. The gnathostome g6pc1 clade was divided into two subclades, each representing the g6pca and g6pcb class genes (see  Supplementary Figure S1 (zs210113_FigS1.pdf)).

Next, an ML tree was constructed based on the nucleotide sequences. In this ML tree, actinopterygian g6pca and g6pcb formed a clade distinct from chondrichthyan and sarcopterygian g6pc1 genes (Fig. 3).

Phylogenetic relationships among amniote g6pc1 genes were inferred from an ML tree constructed based on the nucleotide sequences of the sarcopterygian genes (Fig. 4). In this sarcopterygian ML tree, the amniote g6pc1 genes formed a single clade, suggesting a monophyletic origin of g6pc1 within this group. The g6pc1 clade of amniotes was divided into two subclades. One contained all the mammalian and avian genes, together with at least one gene from each reptilian species. The other subclade consisted exclusively of reptilian genes. We designated the former subclade as “g6pc1-2 class” whereas the latter as “g6pc1-1 class”.

Fig. 2.

Comparative syntenic analyses of g6pc1 genes among sarcopterygian lineages. The synteny around g6pc1 in the sarcopterygian species was compared. While the coelacanth, human, and avian species have only one g6pc1, the western clawed frog and reptiles have two or three paralogs. Regardless of the number of paralogs, g6pc1 genes were flanked by aoc and aarsd1 genes in all species examined.

Fig. 3.

Maximum likelihood phylogenetic tree based on nucleotide sequences of vertebrate g6pc1, g6pc2, and g6pc3. A molecular phylogenetic tree of vertebrate g6pc family genes was constructed using MEGA software (version 7.0.21). The tree was rooted in g6pc3 genes. The scale bar beneath the tree indicates the estimated evolutionary distance units (these settings also apply to Fig. 4). Actinopterygian g6pca and g6pcb formed a clade distinct from chondrichthyan and sarcopterygian g6pc1 genes. This implies that the g6pca and g6pcb classes emerged only in the actinopterygian lineage.

Fig. 4.

Maximum likelihood phylogenetic tree based on nucleotide sequences of sarcopterygian g6pc1, g6pc2, and g6pc3. The nucleotide-based sarcopterygian tree suggested the existence of two g6pc1 classes among amniotes. These were designated g6pc1-1 and g6pc1-2. The former was further divided into two subclasses: g6pc1-1-1 and g6pc1-1-2.

The g6pc1-1 class was further divided into two subclasses, designated as g6pc1-1-1 and g6pc1-1-2. The former contained tuatara, squamates, and one paralog from each of the chelonian species. The latter contained crocodilian paralogs and one paralog each from the green sea turtle and the western painted turtle.

Comparison of genomic environment upstream of g6pc1

Genomic regions upstream of the transcription start site were compared to infer phylogenetic relationships among the reptilian genes (Fig. 5).

When the painted turtle g6pc1-1-1 was used as a query, a similar region was detected in sea turtle g6pc1-1-1 and softshell turtle g6pc1-1 (Fig. 5A). On the other hand, regions similar to painted turtle g6pc1-1-2 were detected in sea turtle g6pc1-1-2 and crocodilian g6pc1-1 (Fig. 5B). In both comparisons, no similarity was detected between the g6pc1-1-1 and g6pc1-1-2 subclasses.

Hepatic gene expression of g6pc1 genes and pck2 in the Japanese gecko

Fasting had no significant effect on the expression of g6pc1 genes. In contrast, pck2 expression was significantly downregulated (Fig. 6A). Comparison of gene expression in free-living female geckos between seasons revealed significantly lower expression of g6pc1 and pck2 in April than in September. In contrast, males showed no significant change between seasons (Fig. 6B). Male geckos administered E2 intraperitoneally showed significantly lower expression levels of g6pc1-1 and pck2 compared to those administered sesame oil alone. In contrast, the expression of g6pc1-2 showed no significant change between the E2 treatment and control groups (Fig. 6C).

Fig. 5.

Comparison of genomic sequences upstream of g6pc1. Nucleotide sequences upstream of the g6pc1 paralogs (–1 kbp from the transcription start site) were aligned and compared using the Vista package. Sequence similarities of over 50% are represented as peaks. Conserved non-coding sequences (CNSs) detected by the program are highlighted. (A) When the painted turtle g6pc1-1-1 was used as a query, similar regions were detected in sea turtle g6pc1-1-1 and softshell turtle g6pc1-1. (B) When the painted turtle g6pc1-1-2 was used as a query, similarity was detected in sea turtle g6pc1-1-2 and crocodilian g6pc1-1. No similarity was detected between g6pc1-1-1 and g6pc1-1-2 subclasses in any of the comparisons described above.

Search for hormone responsive elements upstream of g6pc1 genes

A single ERE (AGGTCA) was present – 928bp upstream of the g6pc1-1 transcription start site. In contrast, no ERE or its half site was present at least 1kbp upstream of g6pc1-2. No PRE was present upstream of either of the g6pc1 genes.

Conservation of functional sites among sarcopterygian g6pc1 paralogs

As mentioned in the Introduction, the primary function of g6pc1 is the dephosphorylation of G6P. In this study, we investigated whether functionally important active center sites and a phosphatase motif (Chou and Mansfield, 2008) are conserved among sarcopterygian g6pc1 genes. For this purpose, the amino acid sequences of g6pc1 from species representing each sarcopterygian lineage were aligned. Alignment revealed that all sarcopterygian genes conserved five active center sites and a phosphatase motif (Fig. 1). Such conservation of sites suggests that all genes retain the catalytic activity of G6P dephosphorylation.

Fig. 6.

Quantitative RT-PCR analyses of Japanese gecko g6pc1 and pck2 in the liver. The results are shown as the mean ± SEM, and the expression levels of the control group (fasting and E2 administration) and males captured in April (seasonal variation) are expressed as 1. * P < 0.05, ** P < 0.01 (Mann–Whitney U test). (A) Gene expression in response to fasting for 3 weeks. The expression of g6pc1 genes showed no significant difference between groups. In contrast, pck2 expression was significantly downregulated by fasting. (B) Seasonal changes in gene expression were examined in males and females. Although no significant seasonal change was observed in males, the expression of both g6pc1 genes and pck2 in females was significantly lower in April than in September. (C) Gene expression in response to E2 treatment. The expression of g6pc1-1 and pck2 was significantly lower in the E2 treatment group, whereas that of g6pc1-2 was not significantly different between groups.

Evolutionary process of duplicated g6pc1 genes

The evolutionary process of g6pc1 was investigated by drawing phylogenetic trees, syntenic analyses, and comparing the upstream sequences. Three molecular phylogenetic trees were constructed in this study, of which two are presented. Although the bootstrap robustness of internal nodes is not always good, insightful information can be obtained from the other analyses.

Gene duplication around actinopterygian-sarcopterygian split

In the amino acid-based NJ tree, the g6pc1 clade was divided into the g6pca and g6pcb subclades (see  Supplementary Figure S1 (zs210113_FigS1.pdf)). This tree topology was in accordance with a previous report (Marandel et al., 2017), which suggested that 2R-WGD in the ancestral gnathostome gave rise to subclades. However, the nucleotide-based ML tree constructed in this study did not support a sister relationship between g6pca and g6pcb. Instead, the actinopterygian g6pca and g6pcb genes formed a distinct clade from the chondrichthyan and sarcopterygian g6pc1 genes (Fig. 3). This tree topology suggested that g6pca and g6pcb emerged in the actinopterygian lineage after the actinopterygian-sarcopterygian split. Therefore, the evolutionary origins of g6pca and g6pcb may need to be reconsidered.

Basal sarcopterygian and amphibian paralogs

In the sarcopterygian ML tree, all g6pc1 genes except that of coelacanth formed a monophyletic clade, which was designated as “tetrapod g6pc1 clade” (Fig. 4). Syntenic analyses revealed only a single g6pc1 gene in the coelacanth, whereas multiple genes were present among the tetrapod species (Fig. 2). This observation suggests that the sarcopterygian g6pc1 genes have emanated from a single ancestral gene.

Among the amphibian species, Nanorana parkeri was found to possess only a single g6pc1 gene. In contrast, other species, including anurans and caecilians, were found to possess two to four genes (Fig. 4), suggesting lineage-specific gene duplication.

Gene duplication in the ancestral amniotes

Our syntenic analyses revealed that mammals possess only a single g6pc1 gene, whereas reptiles possess two to three genes (Fig. 2). This distribution of genes within amniotes suggests that gene duplication occurred in the ancestral reptilian lineage after the mammalian-reptilian split. Molecular phylogeny supports this hypothesis. In both ML trees drawn in this study, the amniote g6pc1 clade bifurcated between the mammalian clade and reptilian-avian clade (Figs. 3, 4).

Gene duplication and loss among reptilian lineages

In the sarcopterygian ML tree, the reptilian-avian clade was divided into two subclades, which were designated as g6pc1-1 and g6pc1-2 classes (Fig. 4). Squamates, crocodilians, and softshell turtles were found to possess two genes, one from each subclass. The tuatara had one more g6pc1-2 class gene, which was likely produced by lineage-specific duplication after divergence from squamates (Fig. 4).

The sea turtle and painted turtle were found to possess three g6pc1 genes, two from the g6pc1-1 class and one from the g6pc1-2 class. The two g6pc1-1 genes were designated g6pc1-1-1 and g6pc1-1-2 (Fig. 4). According to Crawford et al. (2015), these two species and the softshell turtle belong to two sister superfamilies: the former two to Durocryptodira, the latter one to Trionychia. Considering the sister-relationship between these superfamilies, the simplest explanation for the origin of the g6pc1-1-1 and g6pc1-1-2 subclasses is the duplication of the ancestral g6pc1-1 gene in the durocryptodiran lineage (Fig. 7A). However, this hypothesis is not supported by molecular phylogeny. In the sarcopterygian ML tree, durocryptodiran g6pc1-1-1 and g6pc1-1-2 were not immediate siblings. Instead, the former was sister to the softshell turtle g6pc1-1-1, whereas the latter was sister to the crocodilian g6pc1-1 (Fig. 4). An alternative hypothesis that can explain such a tree topology is as follows: First, g6pc1-1 duplicated in the common ancestor of chelonians and crocodilians, giving rise to g6pc1-1-1 and g6pc1-1-2 subclasses. While durocryptodirans retained both subclasses, the softshell turtle and crocodilians lost one of the two subclasses (Fig. 7B). Support for this alternative hypothesis was obtained by comparing the genomic environment. Regions similar to durocryptodiran g6pc1-1-1 were detected in softshell turtle g6pc1-1 but not in durocryptodiran g6pc1-1-2 (Fig. 5A). In contrast, regions similar to g6pc1-1-2 were detected in crocodilian g6pc1-1, but not in durocryptodiran g6pc1-1-1 (Fig. 5B). These results suggest a close relationship between g6pc1-1-1 and softshell turtle g6pc1-1 and between g6pc1-1-2 and crocodilian g6pc1-1.

Loss of g6pc1-1 class from avian lineage

Syntenic analyses revealed only a single g6pc1 gene in the avian species (Fig. 2). This gene belonged to the g6pc1-2 class in the sarcopterygian ML tree (Fig. 4). Although the zebra finch was previously reported to possess two g6pc1 genes (Marandel et al., 2017), only a single gene was present in the current version of the genome assembly (bTaeGut2.pat.W.v2). The presence of two g6pc1 classes in crocodilians suggests the loss of g6pc1-1 from the avian lineage after crocodilian-avian split. Interestingly, mammals possess only one g6pc1 gene. Since mammals and aves are two amniote lineages that have acquired endothermy independently from each other, acquisition and loss of the duplicated g6pc1 gene could be related to thermoregulatory strategies (Fig. 8).

Fig. 7.

Hypotheses regarding the evolution of reptilian g6pc1. (A) Presence of two g6pc1-1 subclasses (g6pc1-1-1 and g6pc1-1-2) only in the durocryptodiran turtles suggests duplication within this lineage. (B) In contrast, the sarcopterygian ML tree and similarity of genomic environment suggest duplication in the common ancestor of crocodilians and chelonians, followed by lineage-specific loss of subclasses.

Transcriptional regulation of gecko g6pc1 genes and implications for physiological significance

In this study, the functional divergence between g6pc1 genes was investigated, focusing on transcriptional regulation in the liver. This is because the transcription of g6pc1 is regulated by multiple hormonal and nutritional stimuli, which in turn results in altered glucose production from gluconeogenic organs (Caseras et al., 2002; Salgado et al., 2004; Hutton and O'Brien, 2009; Yabaluri and Bashyam, 2010; Rui, 2014). In addition, paralog-specific responses to nutritional status have been reported in rainbow trout (Marandel et al., 2015). We also investigated the gene expression of phosphoenolpyruvate carboxykinase 2 (pck2) because pck2 catalyzes the first irreversible step of gluconeogenesis and thus is expected to reflect hepatic gluconeogenic activities.

In the current study, neither of the g6pc1 genes showed significant changes in expression after 3 weeks of fasting. In contrast, pck2 expression was significantly downregulated after fasting (Fig. 6A). Downregulation of pck under prolonged fasting has also been reported in other vertebrates, such as the large yellow croaker (Larimichthys crocea) (Qian et al., 2016). In mammals and birds, gluconeogenesis is repressed within a day or two because excessive breakdown of proteins is harmful (Berg et al., 2002; Wang et al., 2006). Such protein sparing could also be the reason for the downregulation of the gecko pck2 gene.

Fig. 8.

Distribution of g6pc1 paralogs and thermoregulatory strategies among amniotes. Among amniotes, mammals and aves acquired endothermy independently of each other. These endothermic lineages possess only a single g6pc1-2 class paralog. In contrast, all ectothermic reptiles retained both the g6pc1-1 and g6pc1-2 class paralogs. Such a distribution of genes suggests that the retainment of more than one g6pc1 is beneficial only to ectotherms.

Stable expression of gecko g6pc1 genes under prolonged fasting might reflect reliance on glycogen for energy storage, because g6pc1 catalyzes the most downstream reaction shared by gluconeogenesis and glycogenolysis. Although glycogen has been considered to serve only as a short-term energy storage in mammals (Wang et al., 2006), it serves as a long-term energy source in other vertebrates, such as side-blotched lizard (Uta stansburiana) (Zani et al., 2012) and common carp (Cyprinus carpio) (Wang et al., 2006).

Seasonal gene expression changes were significant only in females, which showed lower expression of both g6pc1 and pck2 in April than in September (Fig. 6B). This implies the involvement of sex-related factors. One possible factor is estrogen. Rise of blood E2 levels during vitellogenesis is reported from two close relatives of the Japanese gecko: the leopard gecko (Rhen et al., 2000) and Madagascar ground gecko (Weiser et al., 2012). Since female Japanese geckos undergo vitellogenesis from April to July on Honshu Island (Ikeuchi, 2004), an effect of estrogen was expected. E2 administration resulted in the downregulation of g6pc1-1 and pck2 expression, but not g6pc1-2 expression (Fig. 6C). To infer the molecular basis behind such paralog-specific responses of g6pc1, the ERE was searched within genomic regions up to 1 kbp upstream of the g6pc1 transcription start site. A single ERE half-site was present –928 bp upstream of g6pc1-1, whereas no ERE was present upstream of g6pc1-2. These results suggest that the paralog-specific E2 response is determined by the distribution of the ERE upstream of each gene.

As g6pc1-2 expression in female geckos was lower in April than in September, seasonal and sex-specific factors other than estrogen were expected to regulate this gene. One possible candidate for this was progesterone. Blood progesterone levels have been reported to increase during the reproductive period in other gecko species, including the leopard gecko and Madagascar ground gecko (Rhen et al., 2000; Weiser et al., 2012). However, a PRE was not present within at least 1 kbp upstream of either of the g6pc1 genes. This result implies that progesterone is not involved in the regulation of g6pc1 genes. Further investigation is required to fully understand the seasonal expression of g6pc1-2.

The advantage of such paralog-specific responses to E2 could be related to the maintenance of fertility. The downregulation of g6pc1-1 and pck2 expression by E2 suggests that this hormone represses gluconeogenesis. Similar observations have been reported for both mammalian and non-mammalian vertebrates (Petersen et al., 1983; Qiu et al., 2017). Such downregulation has been proposed to be beneficial for the maintenance of glucose homeostasis under normal conditions (Qiu et al., 2017). However, this can adversely affect female reproduction. Using medaka (Oryzias latipes) as a model, Hasebe et al. (2017) revealed that the firing activity of female GnRH1, which is essential for reproduction, is induced only when blood glucose levels are sufficiently high. That report implies that maintenance of blood glucose levels is essential for reproduction. However, an increase in blood estrogen levels during vitellogenesis can downregulate gluconeogenesis, resulting in lower blood glucose levels and impaired fertility. Under such conditions, estrogen-insensitive g6pc1-2 can maintain glucose production and safeguard female fertility. One important feature of g6pc1 is that it facilitates not only gluconeogenesis but also glycogenolysis. While other gluconeogenic enzymes, such as pck2, only regulate reactions in the gluconeogenesis pathway, g6pc1 catalyzes the common and last reaction shared by gluconeogenesis and glycogenolysis. This glycogenolytic function enables g6pc1 to safeguard glucose production even when gluconeogenesis is downregulated.

The “safeguarding” role of duplicated g6pc1 genes could also explain why non-avian reptiles retained more than one g6pc1 gene (Fig. 8). According to Wieser (1985), ectotherms invest a larger portion of metabolizable energy into reproduction than endotherms do, and are susceptible to mortality from excess investment in reproduction. As a result, ectotherms need to take great care when partitioning energy between their own survivorship and reproduction. Duplicated g6pc1 genes may aid in energy partitioning by fine-tuning glucose production. As described above, g6pc1 catalyzes the last reaction shared by gluconeogenesis and glycogenolysis. Thus, this gene is in a position to singlehandedly affect the net glucose production. By retaining more than one g6pc1 gene and putting them under independent control, ectotherms, in theory, can enjoy greater flexibility in regulating glucose production than endotherms, which only possess a single g6pc1.