Phosphoenolpyruvate Carboxykinase 1 Gene (Pck1) Displays Parallel Evolution between Old World and New World Fruit Bats

Bats are an ideal mammalian group for exploring adaptations to fasting due to their large variety of diets and because fasting is a regular part of their life cycle. Mammals fed on a carbohydrate-rich diet experience a rapid decrease in blood glucose levels during a fast, thus, the development of mechanisms to resist the consequences of regular fasts, experienced on a daily basis, must have been crucial in the evolution of frugivorous bats. Phosphoenolpyruvate carboxykinase 1 (PEPCK1, encoded by the Pck1 gene) is the rate-limiting enzyme in gluconeogenesis and is largely responsible for the maintenance of glucose homeostasis during fasting in fruit-eating bats. To test whether Pck1 has experienced adaptive evolution in frugivorous bats, we obtained Pck1 coding sequence from 20 species of bats, including five Old World fruit bats (OWFBs) (Pteropodidae) and two New World fruit bats (NWFBs) (Phyllostomidae). Our molecular evolutionary analyses of these sequences revealed that Pck1 was under purifying selection in both Old World and New World fruit bats with no evidence of positive selection detected in either ancestral branch leading to fruit bats. Interestingly, however, six specific amino acid substitutions were detected on the ancestral lineage of OWFBs. In addition, we found considerable evidence for parallel evolution, at the amino acid level, between the PEPCK1 sequences of Old World fruit bats and New World fruit bats. Test for parallel evolution showed that four parallel substitutions (Q276R, R503H, I558V and Q593R) were driven by natural selection. Our study provides evidence that Pck1 underwent parallel evolution between Old World and New World fruit bats, two lineages of mammals that feed on a carbohydrate-rich diet and experience regular periods of fasting as part of their life cycle.


Introduction
Food limitation is a common physiological challenge faced by animals and leads to a risk of death due to starvation. Prolonged fasting, and starvation, causes a reduction in blood glucose We speculated that Pck1 experienced adaptive evolution in OWFBs and NWFBs in association with their acquired high carbohydrate diets. To test our hypothesis, we obtained the coding region of Pck1 from 20 species of bats (including five OWFBs species and two NWFBs species), and assessed whether these sequences display molecular evolutionary patterns consistent with adaptive evolution. In addition, we also examined whether the Pck1 gene has experienced convergent and/or parallel evolution between the OWFBs and NWFBs at the amino acid level. Scotophilus heathi was captured in Jul. 2011 from a building (23°1 0 N, 114°56 0 E) located in Guangdong province, China. Artibeus lituratus, Leptonycteris yerbabuenae and Pteronotus parnellii were captured in Apr. 2010 in Mexico. Bats captured in these locations were permitted by the local Protection and Research Center. No bat species used for this study are considered endangered. Bats were sacrificed by decapitation, quickly after capture, and all efforts were made to minimize potential pain and suffering.
To expand our dataset across mammals we incorporated published Pck1 sequences from nine other mammalian species from NCBI and Ensembl: Homo sapiens (NM_002591), Mus musculus (NM_011044), Rattus norvegicus (NM_198780), Equus caballus (XM_001489771), Bos taurus (NM_174737), Sus scrofa (NM_001123158), Canis familiaris (NM_001197143), Loxodonta Africana (XM_003419943), and Monodelphis domestica (XM_001377770). Details on all species and their corresponding accession numbers are listed in S1 Table. Isolation, Amplification and Sequencing of Pck1 Coding Sequences For the 15 bat species used to isolate new sequences, total RNA was isolated from liver tissue (stored at -80°C) using the Trizol reagent. Following the standard protocol, 5 μg of total RNA was reverse transcribed into cDNA using the Reverse Transcriptase kit. Primers (F: ATGCCTCCTCAGCTGCAAAACG and R: TACATCTGGCTGACTCTCTGCCTC) were designed to amplify the coding sequences of Pck1. All PCR products were separated using 1% agarose gels, with the correct sized fragments purified with Gel Extraction Kits, ligated into pGEM-T easy vector, cloned, and sequenced using Terminator kits on an ABI 3730 DNA sequencer. To avoid artifacts, multiple clones (>5) were sequenced for each species.

Phylogenetic Reconstruction
Pck1 nucleotide sequences from 29 mammalian species were aligned using ClustalX [33] and checked for accuracy by eye, with the coding sequences translated into amino acids using MEGA5 [34]. A Bayesian phylogenetic tree was reconstructed based on the aligned nucleotide sequences using MrBayes 3.1.2 [35] with the GTR+I+Γ nucleotide substitution model selected as the best-fitting model by jModelTest0.1 [36]. For the Bayesian analysis, we performed 10,000,000 generations of MCMC and sampled every 100 generations, with the first 2,000,000 generations discarded as burn-in. All other options and priors were the default settings of the MrBayes 3.1.2 software. The standard deviations of split frequencies were stably below 0.01 after 2,000,000 generations of MCMC performances.

Molecular evolution analyses
Before beginning the molecular evolution analyses, we first tested the sequences for evidence of recombination, since recombination may lead to the incorrect results in tests for positive selection. We used GARD [37] from the HyPhy package [38] to detect evidence for recombination breakpoints and tested their statistical significance. For tests of selection we used a topology for the phylogenetic tree that was based on the accepted species relationships [39][40][41][42]. For positive selection in Pck1, we estimated the rates of nonsynonymous (d N ) and synonymous (d S ) substitutions using CODEML from the PAML package [43]. We first used the two-ratio model, where the d N /d S ratio (termed omega or ω) was allowed to be different on the foreground and background branches. To detect differences in the selective pressure acting on the ancestral branches leading to OWFBs and NWFBs, each of these two groups were separately set as the foreground branch. The one ratio-model, where the omega values were equal among branches, served as the null hypothesis for these comparisons. We also employed the free-ratio model, where the value of d N /d S varies among all branches, to explore the overall selection pressure acting on Pck1 in the 29 mammalian species used in this study [44]. Similarly, the one-ratio model served as the null expectation for this comparison. "PairwiseRelativeRate.bf" from the HyPhy package [38] was also used to detect whether there were significant differences in the rates of synonymous and nonsynonymous substitutions between the frugivorous and insectivorous bats (see S1 Text for details).
The branch-site model A (test1 and test2) was then applied to detect positively selected changes along particular branches [45]. In this model, the branch of interest is set as the foreground, with the remaining branches set as the background. We conducted both test1 and test2 of the branch-site models [46] on the ancestral branches for the OWFBs and NWFBs. In both tests, the branch-site model A is the alternative hypothesis. In test1, the null hypothesis is the M1a (Nearly Neutral) model. For test2, the null hypothesis is the branch-site null model. Results from the alternative and null hypothesis were then compared in likelihood ratio tests (LRTs).
Ancestral sequences were reconstructed using CODEML from the PAML package [47], and the amino acid substitutions on each branch were then inferred. The identity and direction of the substitutions was also assessed using the maximum parsimony method with Mesquite 2.74 [48]. To further explore the heterogeneous selection pressure acting across codons of the Pck1 gene in frugivorous bats, a sliding window analysis of the ω values was performed using the program SWAAP 1.0.2 [49] with the Nei and Gojobori method [50]. Window and step sizes for this analysis were set as 90 and 36 nucleotides, respectively.
To detect evidence of parallel or convergent amino acid substitutions between OWFBs and NWFBs, we used a method described by Castoe et al. implemented in the codeMLancestral package [51] (see S1 Text). To test whether the convergence detected between pairs of focal branches was significant, we used a method described by Rossiter et al [52] (see S1 Text). We also used CONVERGE2 [53] to test whether the number of observed parallel or convergent amino acid substitutions shared by the OWFBs and NWFBs branches were significantly different from those expected by random chance. Finally, we examined whether the locations of the positively selected, convergently evolved and other specific sites in the protein sequences were associated with particular functional domains. The protein domains were defined from previous research [54,55].

Results
Our Pck1 sequence dataset contained 29 taxa, including five Old World fruit bats (Pteropodidae), two New World fruit bats (Phyllostomidae) and 13 insectivorous bats. The alignment of the Pck1 coding sequence spanned 1818 nucleotides (representing~97% of the complete coding sequence) and equates to 606 amino acids. No insertions, deletions or stop codons were identified in the isolated coding sequences from the 20 bat species, and no evidence for recombination breakpoints were detected in our dataset.
To detect evidence for changes in the selection pressure acting on Pck1 in OWFBs and NWFBs, we set either the ancestral branch to OWFBs or NWFBs as the foreground branch and performed two-ratio model tests on these branches. Our results showed that the two-ratio model with the ancestral branch to OWFBs as the foreground was a better fit than the oneratio model (null hypothesis) [likelihood ratio test (LRT) statistic (2Δℓ) = 3.99, df = 1, P = 0.046] (Table 1), and had an omega value (d N /d S ) on the ancestral branch leading to OWFBs that was higher than on the background (0.1090 vs. 0.0676, respectively). These results indicate that a change in selection pressure on the Pck1 gene occurred on the OWFB lineage. However, the two-ratio model test with the ancestral branch for NWFBs as the foreground was not significantly different from the one-ratio model (2Δℓ = 3.09, df = 1, P = 0.079) ( Table 1), which implied that no change in selection pressure occurred on this NWFB lineage.
Pairwise relative rate tests were also conducted to determine whether there were significant differences in the rates of synonymous and nonsynonymous substitutions between frugivorous and insectivorous bats. Our results showed that a significant difference in the rate of  nonsynonymous substitutions occurs between OWFBs and insectivorous bats; however, the same result was not found for synonymous substitutions. These results from the pairwise relative rate tests indicated that a change in selection pressure on the Pck1 gene occurred in OWFBs, and that this change mainly effected the rate of nonsynonymous substitutions ( Table 2). Branch-site model A test was performed to detect evidence for positive selection on the ancestral branches of OWFBs and NWFBs. The result of branch-site model A revealed 13 sites on the ancestral branch of OWFBs with elevated omega values (10G, 13G, 62A, 66L, 95V, 188D, 270V, 349E, 352R, 375E, 534V, 578L and 591V, amino acid positions refer to the human PEPCK1 sequence), and the fit of branch-site model A was better than the M1a model (2Δℓ = 7.95, df = 2, P = 0.0188; test1) (Table 3). However, the in-depth test for positive selection (test2) did not reveal a similar result, as the branch-site model A of OWFBs was not significantly better than the null branch-site null model (2Δℓ = 0.0, df = 1, P = 1.0; test2) ( Table 3). Since test1 cannot distinguish relaxed constraints from positive selection [46], we were unable to rule out the possibility of relaxed evolution. No evidence for a change in selection pressure was found in the ancestral branch of NWFBs with both test1 and test2, as the alternative hypothesis of both the M1a model and branch-site null model were not significantly better than the null models.
To identify amino acid changes on particular branches, we use the software Mesquite 2.74 [48] to reconstruct the ancestral amino acid sequences. These reconstructions showed that there were six specific amino acid replacements (A13G, S252G, K349E, F578L, D591V, and I608V) on the ancestral branch of OWFBs (Fig. 2), and four (A13G, K349E, F578L, and D591V) of them were detected by branch-site model A as having elevated omega values (Table 3). To further explore the selection pressure changes among Pck1 codons, we calculated the values of ω among OWFBs, NWFBs, and insectivorous bats using SWAAP 1.0.2. The results of this analysis showed that the estimated ω values for OWFBs are much larger than those of insectivorous bats, with a maximal ω value of 1.9 observed on the nucleotide binding domain 1 (Fig. 3). The ω values of the NWFBs were also larger than those of the insect-eating bats, but not as large as those of the OWFBs.
To examine parallel and/or convergent amino acid substitutions between the OWFB and NWFB ancestral branches, we reconstructed ancestral sequences using PAML [43], with the amino acid replacements on the major branches of the fruit-eating bats inferred. Five parallel amino acid substitutions were shared between OWFBs and NWFBs (N127S, Q276R, R503H, I558V, and Q593R). The total posterior probabilities of divergence versus convergence for focus branch-pairs were shown in Fig. 4. We used two methods to test the probability of acquiring the observed parallel substitutions in the frugivorous branches by chance. Both methods suggest that the probabilities of the substitutions Q276R and I558V, shared between the ancestral branches of OWFBs and A. lituratus, the Q593R substitution, shared between the ancestral branch of NWFBs and C. sphinx, and the R503H substitution, shared between the ancestral branch of the genera Pteropus and L. yerbabuenae, all deviated significantly from the random expectation (Table 4 and Table 5). The N127S substitution, shared between the ancestral branches of OWFBs and L. yerbabuenae, however, was not significant (P > 0.05), suggesting that this amino acid substitution likely occurred in parallel by chance (Table 4 and  Table 5).

Discussion
In this study, we examined the molecular evolution of the Pck1 gene, which plays an important role in the maintenance of glucose homeostasis during fasting in frugivirous bats, which eat carbohydrate-rich diets. As fasting is a regular part of the life cycle of fruit eating bats, they require an enhanced ability to resist hypoglycemia [24,25]. By comparing Pck1 sequences among frugivorous and insectivorous bats, and other mammals, we provide the first data on the molecular evolution of this gene in frugivorous bats. We found that the omega value (d N /d S ) on the ancestral branch of OWFBs is higher than for the remaining branches of the species tree, and that the rate of nonsynonymous substitutions in the Pck1 gene of OWFBs is significantly higher than that of insectivorous bats. Although there was a change in the selection pressure acting on the Pck1 gene on the ancestral branch of NWFBs, the omega value is still quite low indicating that Pck1 is still under purifying selection.
Although the results of test1 of the branch-site model A test revealed 13 amino acid replacement with elevated omega value, we are unable to rule out the possibility that these were due to a relaxation of the selection acting on the pck1 gene, as the result of the explicit test of positive selection (test2) was negative. Test1, by itself, is unable to distinguish between relaxation of constraints and positive selection [46]. However, the likelihood that Pck1 experienced relaxed selection in Old World fruit bats is probably low as PEPCK1 has a very important function across mammals [56]. Even between bacteria and mammals, a number of residues at the active site are conserved and the reaction mechanism for this enzyme in both bacteria and mammals is thought to be same [57]. From our sequenced Pck1 genes we found that the crucial residues in the active site loops, including the R-loop, P-loop and O-loop, are completely conserved.   These results indicated that PEPCK1 retains an important role in glucose homeostasis in frugivore bats [58]. Previous studies have shown that knockout of Pck1 in mice leads to hypoglycemia and death [14], supporting the necessity of PEPCK for life in mammals, and making it unlikely that this gene in OWFBs underwent relaxed evolution. All these results indicate that Pck1 has experienced strong purifying selection not only in fruit bats, but also in all other mammalian species examined. Although the results are congruent with the fact that the PEPCK1protein is very critical for hepatic glycometabolism, there are two possible explanations that may explain the failure to detect positive selection. First, OWFBs diverged from the other Yinpterochiroptera lineages nearly 60 million years ago [40], and previous studies have shown that signals of positive selection can be swamped by a history of purifying selection [46]. Second, the low level of sequence divergence among the sequences limits the power of statistical methods to detect signals of positive selection [59][60][61]. In addition, we only had two NWFB species, which further limits our ability to detect positive selection in this group.
Results from the ancestral reconstruction of the amino acid sequences using the maximum parsimony method with Mesquite 2.74 [48] identified six specific amino acid replacement that occur on the OWFBs ancestral lineage. A change from the basic amino acid lysine to the acidic amino acid glutamic acid (K349E) is located in the PEPCK-specific active site domain (residues 326-425) [54]. The A13G and S252G substitutions occur in the N-terminal domain (residues 1-259) associated with degradation and acetylation of PEPCK1 [62]. F578L, D591V and I608V substitutions lie near the nucleotide-binding domain 2 (residues 426-622) (Fig. 5) responsible for GTP binding, which is necessary for enzymatic activity [55]. Sliding window analysis also found a significantly elevated omega value for the sequences encoding the nucleotide-binding  Fig. 1(B The p value of each test was multiplied by four (branch pairs) to correct for multiple testing. domain 1 (residues 260-325). The elevated omega value is likely due to the I270V substitution, which was detected by the branch-site model A. Replacement of I270V is a OWFBs-specific change among species of the order Chiroptera, and is only shared with the pig (data not shown).
In contrast to OWFBs, we did not find any specific amino acid replacements in the NWFBs. It is interesting to note that other recent molecular evolution studies focusing on diet-related genes in bats also found similar patterns [63][64][65]. For example, GLUT4 which play an important role in glucose homeostasis had two specific amino acid substitutions only in the OWFBs but not in NWFBs [64]. This evolutionary discrepancy in patterns of sequence evolution between OWFBs and NWFBs for the Pck1 gene may have two explanations. First, OWFBs evolved their carbohydrate-rich diet earlier than the NWFBs (28 mya and 20 mya, respectively) [40], thus have had a greater amount of time to accumulated accelerated changes. Second, the specific sequence changes may happen in the regulatory regions in NWFBs, instead of the coding region of these genes, as the promoter can regulate the expression of PEPCK1 [66], change may only occur at the protein abundance and not the protein sequence level.
Although no evidence for positive selection was found in the Pck1 gene of OWFBs and NWFBs, this does not mean that this gene did not undergo adaptive evolution. Parallel evolution between OWFBs and NWFBs is also a strong indicator of adaptive evolution [53]. Previous molecular evolution studies in bats have also found significant evidence for parallel or convergent evolution with no evidence for positive selection [67,68]. A possible explanation for the parallel evolution at the amino acid level between Pck1 genes in OWFBs and NWFBs is molecular adaption to the regular fasting by these fruit bats. It is widely accepted that the hepatic gluconeogenesis plays a critical role in the maintenance of glucose homeostasis during fasting [5,8], and it has been shown that hepatic gluconeogenesis is the primary metabolic mechanism for fasting resistance in the Old World Fruit Bats [31]. PEPCK1, the key enzyme of gluconeogenesis, has been shown in a previous study to have a crucial role in fasting in New World fruit bats [32]. Thus, it is reasonable to assume that the parallel evolution of Pck1 in OWFBs and NWFBs might relate to the evolution of adaptive mechanisms for fasting resistance in fruit bats. Here we found four credible parallel evolution sites (Q276R, R503H, I558V and Q593R) in PEPCK1. All of these sites are located near the nucleotide binding domain, which might influence the binding of GTP to PEPCK1 and thus enzymatic activity [54] (Fig. 5). It is interesting to note that, besides the four convergently evolving sites, the Old World fruit bats also have six other specific amino acid substitutions, which may reflect the fact that the evolutionary history of OWFBs is longer than that of the NWFBs [40].
In conclusion, our results show that the Pck1 gene, which plays an important role in the maintenance of glucose homeostasis during fasting, has undergone parallel evolution between the OWFBs and NWFBs, which provides evidence for adaptive evolution in both lineages, although our evidence for NWFBs is limited to only two species of Phyllostomid bats. Greater sampling of NWFBs is needed to confirm this finding. Apart from gluconeogenesis, PEPCK1 is also the rate-limiting enzyme in glyceroneogenesis, which is responsible for recycling free fatty acid back to triglyceride [69]. Chronic release of free fatty acid by adipose tissue is a critical factor in the development of Type 2 diabetes [70], therefore, PEPCK1 could prevent diabetes by controlling the rate of glyceroneogenesis [70]. Here we found evidence for parallel evolution of the Pck1 gene between OWFBs and NWFBs, which could possibly be explained by its function in glyceroneogenesis. However, additional in-depth studies are needed to address this issue. We speculate that additional deeper studies will reveal that numerous other genes related to carbohydrate metabolism will be found to have experienced adaptive evolution in OWFBs and NWFBs.  Table. Information on species examined in the study. (DOC) S1 Text. Details on some methods used in this study. (DOC)