
==== Front
Genetics
Genetics
genetics
Genetics
0016-6731
1943-2631
Oxford University Press US

37579192
10.1093/genetics/iyad152
iyad152
Investigation
Genome and Systems Biology
AcademicSubjects/SCI01180
AcademicSubjects/SCI01140
Biochemical pathways represented by Gene Ontology-Causal Activity Models identify distinct phenotypes resulting from mutations in pathways
Hill David P The Jackson Laboratory, Bar Harbor, ME 04609, USA

Drabkin Harold J The Jackson Laboratory, Bar Harbor, ME 04609, USA

Smith Cynthia L The Jackson Laboratory, Bar Harbor, ME 04609, USA

Van Auken Kimberly M Division of Biology and Biological Engineering, California Institute of Technology, Pasadena, CA 91125, USA

D’Eustachio Peter Department of Biochemistry and Molecular Pharmacology, New York University Grossman School of Medicine, New York, NY 10016, USA

Baryshnikova A Editor
Corresponding author: Department of Biochemistry and Molecular Pharmacology, New York University Grossman School of Medicine, New York, NY 10016, USA. Email: deustp01@med.nyu.edu
Conflicts of interest The author(s) declare no conflict of interest.

10 2023
14 8 2023
14 8 2023
225 2 iyad15213 7 2023
02 8 2023
29 8 2023
© The Author(s) 2023. Published by Oxford University Press on behalf of The Genetics Society of America.
2023
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Gene inactivation can affect the process(es) in which that gene acts and causally downstream ones, yielding diverse mutant phenotypes. Identifying the genetic pathways resulting in a given phenotype helps us understand how individual genes interact in a functional network. Computable representations of biological pathways include detailed process descriptions in the Reactome Knowledgebase and causal activity flows between molecular functions in Gene Ontology-Causal Activity Models (GO-CAMs). A computational process has been developed to convert Reactome pathways to GO-CAMs. Laboratory mice are widely used models of normal and pathological human processes. We have converted human Reactome GO-CAMs to orthologous mouse GO-CAMs, as a resource to transfer pathway knowledge between humans and model organisms. These mouse GO-CAMs allowed us to define sets of genes that function in a causally connected way. To demonstrate that individual variant genes from connected pathways result in similar but distinguishable phenotypes, we used the genes in our pathway models to cross-query mouse phenotype annotations in the Mouse Genome Database (MGD). Using GO-CAM representations of 2 related but distinct pathways, gluconeogenesis and glycolysis, we show that individual causal paths in gene networks give rise to discrete phenotypic outcomes resulting from perturbations of glycolytic and gluconeogenic genes. The accurate and detailed descriptions of gene interactions recovered in this analysis of well-studied processes suggest that this strategy can be applied to less well-understood processes in less well-studied model systems to predict phenotypic outcomes of novel gene variants and to identify potential gene targets in altered processes.

Genes act in interconnected biological pathways, so single mutations can yield diverse phenotypes. To use the large body of mouse functional gene annotations, we converted human Gene Ontology-Causal Activity Models (GO-CAMs) of glucose metabolism to orthologous mouse GO-CAMs. We then queried phenotypes for mouse genes in these GO-CAMs and identified gene networks associated with discrete phenotypic outcomes due to perturbations of glycolysis and gluconeogenesis. This strategy can be extended to less well-understood processes and model systems to predict phenotypic outcomes.

Reactome Knowledgebase
Gene Ontology-Causal Activity Model (GO-CAM)
Mouse Genome Informatics
glucose metabolism
phenotypes
transcriptional regulation
NIH 10.13039/100000002 U41 HG002273 U24 HG012212 U24 HG012198 U24 HG002223 U41 HG000330 U24 HG011851
==== Body
pmcIntroduction

Gene products mediate networks of causally connected functions to give rise to biological outcomes, revealed as observable phenotypes. Phenotypes correlate with gene functions, and GO biological processes can be predicted by analyzing phenotypic outcomes (Ascensao et al. 2014). Causally connected gene functions can be grouped into biological pathways that represent their coordinated action to achieve a specific outcome. However, phenotypes, even ones conventionally attributed to a single mutated gene, can be the result of interactions of multiple genes. Phenotypic variability among affected individuals described as incomplete penetrance and variable expressivity can be the result of interactions between the mutated gene and other genes with which it normally interacts (Kingdom and Wright 2022). Genes can also display epistatic interactions in which the phenotype caused by a mutation in 1 gene can mask a phenotype of another variant gene in the same pathway (de Visser et al. 2011). Recent studies in large human populations of the range of phenotypes observed in individuals heterozygous for CFTR (cystic fibrosis transmembrane conductance regulator) mutations (Barton et al. 2022), mutations in genes associated with developmental delay (Kingdom et al. 2022), and mutations in genes associated with maturity onset diabetes of the young (Mirshahi et al. 2022) provide striking illustrations of the complex interactions that occur between a mutation in 1 gene and normal genomic polymorphisms resulting in an affected individual's phenotype.

Pathway databases, which model biological processes as networks of reactions, provide useful tools for exploring such interactions. An individual reaction converts input physical entities into modified or relocated output entities. The conversion is mediated by regulatory and catalytic activities of still other entities. A reaction is thereby temporally and causally connected to other reactions that generate its inputs, consume its outputs, or modulate the activities of its catalysts and regulators, and loss of a single protein's function can result in multiple phenotypic consequences determined by the reaction network structure. Understanding this structure thus provides a framework for interpreting phenotypes associated with mutations in genes present in the network.

How generally can this network model be applied to predict phenotypes resulting from the disruption of individual genes given diverse genetic backgrounds? Can we use phenotypes to tease apart subpathways in a network and identify the ones whose disruption gives rise to a discrete phenotype? Here, we used a network of pathways representing glucose and pyruvate metabolism to address these questions. Under normal physiological conditions in vertebrates glucose levels are tightly controlled despite variation in demands for energy and biosynthetic intermediates that depend on glucose metabolism, and variation in glucose availability from dietary intake, glycogen breakdown and de novo glucose synthesis.

The complementary, interconnected processes of glycolysis and gluconeogenesis, in which glucose is catabolized and synthesized, respectively, are central to glucose homeostasis. Glycolysis converts glucose and other monosaccharides to pyruvate and generates ATP and reduced NAD (Prochownik and Wang 2021). Gluconeogenesis consumes pyruvate, NADH, and ATP to generate glucose (Chourpiliadis and Mohiuddin 2022). Lactate (Cori 1981) and amino acids (Felig 1973) are major sources of carbon skeletons that are converted to pyruvate, and there are 3 routes from pyruvate to a common phosphoenolpyruvate intermediate (Fig. 1a). From phosphoenolpyruvate to glucose is a 10-reaction process. Seven reactions are shared by gluconeogenesis and glycolysis: These reactions are reversible under physiological conditions, mediated by the same enzymes, and directed by mass action. The remaining reactions of each process are irreversible under physiological conditions and mediated by different enzymes.

Fig. 1. Reactome pathways. a) Three alternative routes from mitochondrial pyruvate to cytosolic oxaloacetate at the start of the gluconeogenesis pathway (https://reactome.org/PathwayBrowser/#/R-HSA-70263), corresponding to GO-CAM models gluconeogenesis 1 (red highlighting), 2 (yellow), and 3 (magenta). b) Aerobic conversion of pyruvate to acetyl-CoA (green highlighting) and anaerobic conversion to lactate (gray) (https://reactome.org/PathwayBrowser/#/R-HSA-70268). In each case, the process description format shows the substrates and products (chemical entities) of each reaction, the identities of the proteins or complexes that mediate it, and causal connections due to shared products/substrates.

Glycolysis occurs in all tissues in the body. For red blood cells, which lack mitochondria, glycolysis is the only source of ATP. In contrast, gluconeogenesis is mainly confined to hepatocytes and kidney proximal tubules. The irreversible steps of glycolysis and gluconeogenesis are normally coordinately regulated to maximize hepatocyte glycolysis when blood glucose levels are high and in skeletal muscle in response to strenuous exercise, and to minimize hepatocyte glycolysis and maximize hepatocyte and kidney proximal tubule gluconeogenesis in response to fasting and physiological stresses. The pyruvate end-product of glycolysis is metabolized aerobically to acetyl-CoA or anaerobically to lactate (Fig. 1b).

We distinguished 3 groups of genes and regulatory processes for our analysis as follows: those unique to glycolysis, those unique to gluconeogenesis, and those shared between glycolysis and gluconeogenesis, allowing us to ask whether variant forms of genes in each set are associated with distinct or shared phenotypes. Shared phenotypes for genes in 1 set that are distinct from phenotypes of another set would allow us to associate an individual biochemical pathway with specific outcomes. Thus, we could identify genes in the pathway that could potentially be manipulated to modify clinical outcomes.

We have previously computationally generated Gene Ontology-Causal Activity Models (GO-CAMs) of the human processes of glycolysis and gluconeogenesis, and of the metabolism of pyruvate, the end-product of glycolysis, to lactate (anaerobic metabolism) or acetyl-CoA (aerobic metabolism) (Good et al. 2021). The GO paradigm annotates individual gene products with terms that describe their molecular functions, association with cellular components, and involvement in multistep biological processes (The Gene Ontology Consortium 2017). GO-CAM models link these annotations of individual gene products into a causal flow. In GO-CAM models of metabolic processes, this causal flow is represented by enzymatic activities enabled by gene products that are connected by shared input, output molecules or regulatory functions of gene products within the causal chain (Thomas et al. 2019). A GO-CAM differs from the Reactome model of the same process in that the latter employs a more discursive process description data model (Le Novère et al. 2009) to describe processes as networks of reactions mediated by gene products that transform input molecules into output molecules (Gillespie et al. 2022). GO-CAM models derived from Reactome pathways can be used to define sets of genes that function in linear causal pathways and can be used to discriminate genes that function in multiple or branched pathways. A specific feature of GO-CAM models generated in the work described here is that each covers a group of molecular functions corresponding to a single GO biological process, while Reactome models are more variable in size and scope, combining, for example, metabolic and regulatory functions, and merging multiple tissue-specific variant forms of a process.

We have now used computationally-derived human GO-CAM models as templates to manually construct GO-CAM models of the corresponding mouse pathways, teasing apart pathway networks to create a causal flow of GO molecular functions enabled by mouse gene products. We used the gene lists from these models as a basis to interrogate MGD phenotype annotations to determine if mutant forms of mouse genes associated with different mouse GO-CAM models were enriched for different phenotypes (Blake et al. 2021). This approach takes advantage of the wealth of data in the laboratory mouse linking genetic variation to phenotypic variation (http://www.informatics.jax.org/greenbook/).

Laboratory mice are a uniquely rich resource of mammalian phenotypes associated with variations in single genes expressed on diverse genetic backgrounds. A large and diverse collection of mouse genetic variants associated with visible phenotypes has been assembled in the past 120 years by systematically searching for visibly abnormal individuals in inbred populations, by inducing random mutations by mutagens, radiation or insertion of gene trap constructs, and by modern genome editing technologies such as targeted mutagenesis in ES cells or CRISPR/Cas9 methodologies (reviewed in Bello et al. 2021). The variant heritable alleles generated in these ways involve DNA sequence changes including point mutations, insertions of heterologous sequence, chromosomal deletions, and rearrangements affecting several genes, and their phenotypes have been assessed to various extents, by individual investigators and in small and large scale mutagenesis projects (Groza et al. 2023). These disparate phenotype data whether generated by large projects, or by expert curation of published literature, are standardized using terms from the Mammalian Phenotype ontology (https://www.informatics.jax.org/vocab/mp_ontology) (Smith et al. 2005) and to relevant experimental strain genetic backgrounds and evidence, and are integrated in the phenotypes section of the Mouse Genome Informatics resource (https://www.informatics.jax.org/phenotypes.shtml) (Bello et al. 2015). MGI genotype–phenotype data include simple genotypes with causative mutations in a single gene as well as more complex genotypes where mutations in multiple genes may contribute to the phenotype. For analyses here, only the simple genotype associations were used (Bello and Eppig 2016).

Lastly, we asked if we could extend our knowledge mechanistically by using enriched phenotypes for gluconeogenesis to identify transcription factors that are associated with the same phenotypes. We then confirmed that some of the transcription factors we identified are involved in regulation of gluconeogenesis while others are involved in pathways that interact with the gluconeogenic pathway. With this inferential approach, we can begin to systematically expand our pathway analysis to include coordinated transcriptional regulation of the pathway. Combining transcriptional regulation with well-studied allosteric regulation of the enzymes of the pathway, we can begin to better understand the causal connections among genes that result in complex abnormal phenotypes.

Materials and methods

Genes used in GO-CAM models

GO-CAM (Thomas et al. 2019) activity flow (Le Novère et al. 2009) models computed from Reactome gluconeogenesis, glycolysis, pyruvate–lactate, and pyruvate–acetyl-CoA pathways (Hill et al. 2016; Good et al. 2021; Gillespie et al. 2022) assigned functions to 60 human gene products (Table 1). These gene products are grouped in Reactome into sets of paralogs that share a function (e.g. the “glucokinase” set consists of GCK, HK1, HK2, HK3, and HKDC1). We identified high-confidence mouse orthologs of these human gene products: mouse gene products with experimental evidence of conserved function from direct assays (Blake et al. 2021; Ringwald et al. 2022). Supporting evidence, such as references and evidence codes, was added to the mouse models using standard protocols for GO annotation (Balakrishnan et al. 2013; The Gene Ontology Consortium 2017). This search yielded 44 mouse gene products, typically 1 or 2 corresponding to each human set. This search failed for 1, Slc25a1 (MGI:1345283): Despite its sequence similarity to human SLC25A1, we found no experimental evidence to support its role in phosphoenolpyruvate transport in mice and only suggestive indirect evidence in related species. One mouse gene product with no human counterpart, Eno1b (MGI:3648653), was identified: This gene arose by retrotransposition after the divergence of mouse and human lineages from their last common ancestor (Bulusu et al. 2017). The 44 mouse gene products identified in this step were used to construct mouse GO-CAM models. To construct a mouse GO-CAM, each molecular function node in the corresponding human GO-CAM was populated with the 1 or 2 mouse gene products known from direct experimental evidence to enable that function.

Table 1. Mouse and human genes used for analyses.

The part of the Reactome gluconeogenesis model that represented alternative pathways from pyruvate to phosphoenolpyruvate was split into separate mouse GO-CAMs. Separating the Reactome pathways that represent multiple branches of a single process into individual GO-CAMS allows us to compare and contrast “subpathways” as in our analysis of the fates of pyruvate described below. We created 4 pathways to represent glycolysis and gluconeogenesis as follows: canonical glycolysis 1 (mouse)—gomodel:5745387b00001516, gluconeogenesis 1 (mouse)—gomodel:61e0e55600001225, gluconeogenesis 2 (mouse)—gomodel:62d0afa500000132, and gluconeogenesis 3 (mouse)—gomodel:62d0afa500000298 (Fig. 1a and Table 1). For analyses of downstream pathways, we created 2 mouse GO-CAM models based on the Reactome “Pyruvate Metabolism” (R-HSA-70268) pathway, 1 representing the interconversion of pyruvate and lactate—gomodel:633b013300001238 and 1 representing the conversion of pyruvate to acetyl-CoA—gomodel:633b013300001469 (Fig. 1b). All GO-CAM models can be browsed at http://model.geneontology.org/##, where ## is replaced by the model ID (e.g. http://model.geneontology.org/633b013300001238).

Genes used for phenotype enrichment analysis

We extended the list of mouse–human gene product pairs used to construct mouse GO-CAM models by searching MGI for all mouse orthologs of human set members supported by any experimental annotations or sequence orthology annotations consistent with conserved function. The resulting list contains 59 human–mouse pairs, excluding human SLC25A1 (no experimentally verified mouse counterpart) and mouse Eno1b (no human counterpart) (Table 1).

The mouse gene products on this larger list were assembled into 5 sets that were derived from the GO-CAM models and represented genes whose functions are causally connected in each pathway. These GO-CAM-derived “causal” gene sets were used for phenotype enrichment analyses: (1) genes unique to glycolysis, (2) genes unique to gluconeogenesis, (3) genes shared between glycolysis and gluconeogenesis, (4) genes unique to glycolysis plus genes that convert pyruvate to lactate, and (5) genes unique to glycolysis plus genes that convert pyruvate to acetyl-CoA. Enrichment analysis was performed on the gene sets using the VisuaL Annotation Display (VLAD) tool (http://proto.informatics.jax.org/prototypes/vlad/) (Richardson and Bult 2015) using default settings and MGI annotation data complete through 2023 January 26. For comparison of sets 1, 2, and 3, we used the comparative analysis functionality of VLAD. This functionality calculates an enrichment score for each set and plots the relative enrichment significance in graphical format overlaying the structure of the Mammalian Phenotype (MP) ontology (Smith and Eppig 2012). Results are also generated in a tabular format. We hypothesized that genes causally connected in the separate branches would share separate phenotypes with the main branch. We analyzed sets 1, 4, and 5 individually. We filtered the results of each analysis using a false discovery rate (FDR) q < 0.05 cutoff. If a phenotype term, e.g. MP:0013663 “increased myeloid cell number”, and a more specific is_a child of that term, e.g. MP:0002640 “reticulocytosis” both fell below the q < 0.05 cut-off, we considered only the child term in our analyses here. Next, we determined which of these significant phenotypes had a greater number of genes associated with glycolysis and either downstream pathway compared with glycolysis alone or with the other downstream pathway. Complete results from all VLAD analyses are shown in Supplementary Tables 1–5.

Analysis of transcription factors that correspond to phenotypes enriched in the gluconeogenesis gene set

To identify transcription factors that are associated with phenotypes enriched in the gluconeogenesis genes, we searched the MGI “Genes and Markers Query Form” on (https://www.informatics.jax.org/marker—2023 February 2 data) for genes annotated to the GO term “DNA binding transcription factor activity” (GO:0003700) and to the 10 most specific enriched MP terms with a FDR q < 0.05 in the VLAD graph that correspond to specific metabolic or developmental processes: “hypoglycemia” (MP:0000189), “increased liver triglyceride level” (MP:0009355), “decreased circulating triglyceride level” (MP:0002644), “increased liver glycogen level” (MP:0010400), “increased kidney glycogen level” (MP:0031003), “decreased granulocyte number” (MP:0000334), “abnormal midbrain development” (MP:0003864), “abnormal gluconeogenesis” (MP:0003383), “increased circulating triglyceride level” (MP:0001552), and “hepatic steatosis” (MP:0002628).

Results

Creation and validation of mouse GO-CAM models

To determine whether we could associate biochemical pathways with specific phenotypes, we first built mouse GO-CAM models for glycolysis, gluconeogenesis, and pyruvate metabolism. Reactome-derived human GO-CAM models served as templates and were populated with orthologous mouse gene products available from MGI (Table 1). We generated a single mouse pathway for glycolysis. For gluconeogenesis, we generated 3 mouse GO-CAM models, 1 for each of the 3 distinct routes to the formation of cytosolic oxaloacetate in the Reactome model of human gluconeogenesis (https://reactome.org/PathwayBrowser/#/R-HSA-70263), (Fig. 1a and Table 1). We generated models for anaerobic metabolism of pyruvate to lactate and aerobic metabolism to acetyl-CoA, representing branches of pyruvate metabolism in the Reactome human model (https://reactome.org/PathwayBrowser/#/R-HSA-70268) (Fig. 1b and Table 1). For each of the genes in the models, we supported our assertions with standard GO annotation evidence (Balakrishnan et al. 2013). A total of 43 assertions were supported by direct experimental evidence from mouse systems, while 16 were supported by sequence orthology to human or rat proteins whose functions are known from experimental data. Genes associated with glycolysis and gluconeogenesis were classified as involved only in gluconeogenesis, only in glycolysis, or shared by both pathways (Fig. 2 and Table 1).

Fig. 2. Causal connections among gene products that mediate glycolysis, gluconeogenesis, and pyruvate metabolism. Activity flow representation of the pathways of glycolysis, gluconeogenesis, and pyruvate catabolism, showing gene products (Table 1) associated with each process. The enzymes or transporters that mediate the generation of each required small molecule in these processes are connected by arrows to the ones that consume it. Red arrows connect the steps of glycolysis; blue arrows the steps of gluconeogenesis, and green and gray arrows the steps of aerobic and anaerobic pyruvate metabolism. Red shading highlights enzymes and transporters unique to glycolysis; blue shading enzymes and transporters unique to gluconeogenesis; and purple shading enzymes shared by the 2. The green and gray boxes, respectively, indicate the reactions that convert pyruvate to acetyl-CoA or lactate.

Related biological pathways are associated with related but distinct sets of specific phenotypes

To ask whether variant forms of genes associated with a pathway were associated with similar phenotypes, we used the MP ontology terms assigned by MGI to the genes in our mouse GO-CAM pathways (Table 1) to perform comparative enrichment analyses using the VisuaL Annotation Display (VLAD) tool (Richardson and Bult 2015). Fifty-three of 59 genes in our pathways were annotated with MP terms in MGI. We analyzed 3 sets of genes, glycolysis only, gluconeogenesis only, and shared between glycolysis and gluconeogenesis. This VLAD enrichment comparative analysis identified 95 phenotypes that were enriched in at least 1 of the gene sets with a FDR (q) < 0.05 (Table S1). The top scoring phenotypes for comparative enrichment, which fall in the domains of homeostasis, embryonic lethality, and hematopoietic system, and their relationships to one another in the MP ontology are shown graphically in Fig. 3. The 3 colors in the bar in each node represent the relative significance values for each gene set. A tabular summary of the relative enrichment data is shown in Table 2, split into 3 sections as follows: gluconeogenesis only, glycolysis only, and shared. The table shows enrichment results for all MP domains that fit 3 criteria as follows: (1) The term has a FDR < 0.05, (2) the term has no children that are also FDR < 0.05, and (3) when a term or a child passes criteria (1) and (2) in more than 1 set, the terms with the smallest FDR are retained. For example, abnormal hematocrit (shared) is not shown in the table because decreased hematocrit (glycolysis) scores better with these criteria.

Fig. 3. Graphical view of VLAD comparative enrichment analysis for mammalian phenotype annotation data from 2023 January 26. Three gene sets were used for the analysis as follows: glycolysis (green), gluconeogenesis (red), and genes shared between glycolysis and gluconeogenesis (blue) (Table 1). Blue arrows indicate a direct is_a relationship between terms, and black arrows indicate that the terms are indirectly related. As indicated in Table 1, 1 glycolysis-only gene (Pfkp MGI:1891833) and 5 shared genes (Pgk1 MGI:97555, Pgam2 MGI:1933118, Eno1b MGI:3648653, Eno2 MGI:95394, and Eno3 MGI:95395) have not been associated with any variant phenotypes and so were excluded from this analysis. Phenotypes in the domains of homeostasis, embryonic lethality, and hematopoietic system are shown. This graph shows the top 25 phenotypes based on enrichment score for clarity. The complete results of this VLAD analysis are tabulated in Table S1.

Table 2. Enrichment results for the most specific phenotype terms in the comparative analysis.

TermID	Term	Q	# of genes annotated to term	# of genes in set	Query set	Symbols	
Shared set							
MP:0002675	Asthenozoospermia	1.00E+00	0	19	Gluconeogenesis-specific		
1.00E+00	0	9	Glycolysisis-specific		
2.65E−02	3	11	Shared	Eno4, Gapdhs, Pgk2	
MP:0009835	Absent sperm annulus	1.00E+00	0	19	Gluconeogenesis-specific		
1.00E+00	0	9	Glycolysisis-specific		
3.47E−02	1	11	Shared	Eno4	
MP:0005642	Decreased mean corpuscular hemoglobin concentration	1.00E+00	0	19	Gluconeogenesis-specific		
1.00E+00	0	9	Glycolysisis-specific		
3.66E−02	2	11	Shared	Gpi1, Tpi1	
MP:0011609	Decreased glyceraldehyde-3-phosphate dehydrogenase (NAD+) (phosphorylating) activity	1.00E+00	0	19	Gluconeogenesis-specific		
1.00E+00	0	9	Glycolysisis-specific		
2.65E−02	1	11	Shared	Gapdh	
MP:0006205	Embryonic lethality between implantation and somite formation	1.00E+00	0	19	Gluconeogenesis-specific		
1.00E+00	0	9	Glycolysisis-specific		
2.65E−02	3	11	Shared	Eno1, Gapdh, Tpi1	
Glycolysis							
MP:0008859	Abnormal hair cycle catagen phase	1.00E+00	0	19	Gluconeogenesis-specific		
4.93E−02	1	9	Glycolysisis-specific	Pfkl	
1.00E+00	0	11	Shared		
MP:0003731	Abnormal retina outer nuclear layer morphology	1.00E+00	0	19	Gluconeogenesis-specific		
4.93E−02	2	9	Glycolysisis-specific	Hkdc1, Pkm	
1.00E+00	0	11	Shared		
MP:0004952	Increased spleen weight	1.00E+00	0	19	Gluconeogenesis-specific		
4.82E−02	2	9	Glycolysisis-specific	Hk1, Pfkm	
8.36E−02	2	11	Shared	Gpi1, Tpi1	
MP:0013404	Decreased circulating lactate level	1.00E+00	0	19	Gluconeogenesis-specific		
4.65E−02	1	9	Glycolysisis-specific	Pfkm	
1.00E+00	0	11	Shared		
MP:0010035	Increased erythrocyte clearance	1.00E+00	0	19	Gluconeogenesis-specific		
4.65E−02	1	9	Glycolysisis-specific	Pklr	
1.00E+00	0	11	Shared		
MP:0001585	Hemolytic anemia	1.00E+00	0	19	Gluconeogenesis-specific		
6.18E−03	2	9	Glycolysisis-specific	Hk1, Pklr	
1.12E−01	1	11	Shared	Tpi1	
MP:0003059	Decreased insulin secretion	1.00E+00	0	19	Gluconeogenesis-specific		
3.76E−02	2	9	Glycolysisis-specific	Gck, Pfkm	
1.00E+00	0	11	Shared		
MP:0031085	Increased erythrocyte osmotic fragility	1.00E+00	0	19	Gluconeogenesis-specific		
3.08E−03	2	9	Glycolysisis-specific	Pfkm, Pklr	
1.00E+00	0	11	Shared		
MP:0011606	Decreased glucokinase activity	1.00E+00	0	19	Gluconeogenesis-specific		
1.84E−02	1	9	Glycolysisis-specific	Gck	
1.00E+00	0	11	Shared		
MP:0005293	Impaired glucose tolerance	1.00E+00	0	19	Gluconeogenesis-specific		
2.78E−02	3	9	Glycolysisis-specific	Gck, Hkdc1, Pfkm	
1.00E+00	0	11	Shared		
MP:0005639	Hemosiderosis	1.00E+00	0	19	Gluconeogenesis-specific		
4.82E−02	1	9	Glycolysisis-specific	Pklr	
1.00E+00	0	11	Shared		
MP:0010375	Increased kidney iron level	1.00E+00	0	19	Gluconeogenesis-specific		
4.93E−02	1	9	Glycolysisis-specific	Hk1	
1.00E+00	0	11	Shared		
MP:0000245	Abnormal erythropoiesis	1.00E+00	0	19	Gluconeogenesis-specific		
3.77E−02	2	9	Glycolysisis-specific	Hk1, Pklr	
2.28E−01	1	11	Shared	Tpi1	
MP:0002814	Hyperchromasia	1.00E+00	0	19	Gluconeogenesis-specific		
4.36E−02	1	9	Glycolysisis-specific	Hk1	
1.00E+00	0	11	Shared		
MP:0002590	Increased mean corpuscular volume	1.00E+00	0	19	Gluconeogenesis-specific		
4.65E−02	2	9	Glycolysisis-specific	Hk1, Pklr	
8.02E−02	2	11	Shared	Gpi1, Tpi1	
MP:0031318	Enhanced skeletal muscle regeneration	1.00E+00	0	19	Gluconeogenesis-specific		
4.65E−02	1	9	Glycolysisis-specific	Pfkm	
1.00E+00	0	11	Shared		
MP:0000208	Decreased hematocrit	1.00E+00	0	19	Gluconeogenesis-specific		
1.11E−02	3	9	Glycolysisis-specific	Hk1, Pfkm, Pklr	
9.30E−02	2	11	Shared	Gpi1, Tpi1	
MP:0002640	Reticulocytosis	1.00E+00	0	19	Gluconeogenesis-specific		
1.94E−03	3	9	Glycolysisis-specific	Hk1, Pfkm, Pklr	
3.47E−02	2	11	Shared	Gpi1, Tpi1	
MP:0005344	Increased circulating bilirubin level	3.16E−01	1	19	Gluconeogenesis-specific	Slc37a1	
4.60E−03	3	9	Glycolysisis-specific	Hk1, Pfkm, Pklr	
5.87E−02	2	11	Shared	Gpi1, Tpi1	
MP:0000874	Irregular external granule cell layer thickness	1.00E+00	0	19	Gluconeogenesis-specific		
1.34E−02	1	9	Glycolysisis-specific	Hk2	
1.00E+00	0	11	Shared		
MP:0006353	Increased glycosylated hemoglobin level	1.00E+00	0	19	Gluconeogenesis-specific		
4.52E−02	1	9	Glycolysisis-specific	Gck	
1.00E+00	0	11	Shared		
Gluconeogenesis						
MP:0000334	Decreased granulocyte number	3.08E−02	3	19	Gluconeogenesis-specific	G6pc3, Got1, Slc37a4	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0003864	Abnormal midbrain development	2.44E−02	2	19	Gluconeogenesis-specific	Got2, Mdh2	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0001731	Abnormal postnatal growth	4.63E−02	5	19	Gluconeogenesis-specific	G6pc, G6pc3, Pck1, Slc25a12, Slc37a4	
4.40E−01	1	9	Glycolysisis-specific	Gck	
1.00E+00	0	11	Shared		
MP:0003383	Abnormal gluconeogenesis	3.08E−02	2	19	Gluconeogenesis-specific	Pcx, Slc25a13	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0013405	Increased circulating lactate level	7.85E−03	2	19	Gluconeogenesis-specific	Pck1, Pcx	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0002575	Increased circulating ketone body level	1.62E−02	2	19	Gluconeogenesis-specific	Pck1, Pcx	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0000189	Hypoglycemia	9.96E−07	6	19	Gluconeogenesis-specific	Fbp1, G6pc, Pck1, Pcx, Slc25a13, Slc37a4	
1.27E−01	1	9	Glycolysisis-specific	Gck	
1.00E+00	0	11	Shared		
MP:0002628	Hepatic steatosis	5.28E−03	4	19	Gluconeogenesis-specific	G6pc, Pck1, Pcx, Slc25a13	
1.97E−01	1	9	Glycolysisis-specific	Gck	
3.26E−01	1	11	Shared	Aldob	
MP:0009355	Increased liver triglyceride level	6.31E−04	4	19	Gluconeogenesis-specific	G6pc, Pck1, Pcx, Slc25a13	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0010400	Increased liver glycogen level	1.61E−02	2	19	Gluconeogenesis-specific	G6pc, Slc37a4	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0002644	Decreased circulating triglyceride level	3.58E−02	3	19	Gluconeogenesis-specific	G6pc2, Pck1, Pck2	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		
MP:0001552	Increased circulating triglyceride level	7.05E−03	4	19	Gluconeogenesis-specific	G6pc, Pck1, Pcx, Slc37a4	
2.09E−01	1	9	Glycolysisis-specific	Gck	
1.00E+00	0	11	Shared		
MP:0031003	Increased kidney glycogen level	5.28E−04	2	19	Gluconeogenesis-specific	G6pc, Slc37a4	
1.00E+00	0	9	Glycolysisis-specific		
1.00E+00	0	11	Shared		

The 3 gene sets were associated with distinct enriched phenotypes. The gluconeogenesis genes had more significant enrichment in categories representing “abnormal triglyceride levels”, “abnormal glycogen levels”, “hypoglycemia”, “ketosis”, and growth and development. The glycolysis genes had more significant enrichment in categories representing glucose homeostatic phenotypes as well as erythrocyte (in)stability and downstream consequences such as iron levels and bilirubin levels. Shared genes were more significantly enriched for sperm and lethal phenotypes.

The higher levels of the MP ontology hierarchy in Fig. 3 show that top scoring abnormal phenotypes for genes unique to glycolysis tend to affect erythropoiesis, while those of gluconeogenesis tend to affect energy metabolism, notably of glycogen and triglycerides.

By convention, glycolysis ends with production of pyruvate. Pyruvate is further metabolized aerobically to form acetyl-CoA or anaerobically to form lactate (Figs. 1b and 2a) (Prochownik and Wang 2021). Given the distinct phenotypes associated with variation in genes associated with glycolysis and gluconeogenesis, we next tested if we could discriminate between phenotypes associated with these 2 fates of pyruvate. We compared 3 gene sets as follows: genes that were unique to glycolysis (i.e. not in gluconeogenesis) and these glycolysis genes in combination with those either for aerobic or anaerobic pyruvate metabolism. We analyzed the 3 gene sets using VLAD and compared phenotypes that were enriched with a FDR < 0.05 (Table 3) by counting the number of genes associated with each phenotype.

Table 3. Comparison of number of genes annotated to phenotypes in downstream pathways.

ID	Terma	Lactate	Glycolysis	Acetyl-CoA	
Phenotypes with more genes in the glucose to lactate pathway (anaerobic) than glucose to pyruvate	
MP:0002447	Abnormal erythrocyte morphology	5 (5.24E−02)	4 (4.55E−02)	4 (1.51E−01)	
MP:0002875	Decreased erythrocyte cell number	3 (3.71E−02)	2 (6.79E−02)	2 (1.51E−01)	
MP:0000208	Decreased hematocrit	4 (4.21E−03)	3 (1.11E−03)	3 (3.37E−02)	
MP:0005563	Abnormal hemoglobin content	4 (8.64E−03)	3 (1.74E−02)	3 (5.74E−02)	
MP:0002874	Decreased hemoglobin content	3 (2.69E−02)	2 (5.21E−02)	2 (1.27E−01)	
MP:0002590	Increased mean corpuscular volume	3 (2.03E−02)	2 (4.65E−02)	2 (1.07E−01)	
MP:0013707	Abnormal hematopoietic precursor cell morphology	3 (3.82E−02)	2 (7.01E−02)	2 (1.51E−01)	
MP:0001585	Hemolytic anemia	3 (4.06E−04)	2 (6.18E−03)	2 (1.74E−02)	
MP:0002640	Reticulocytosis	4 (2.10E−04)	3 (1.94E−03)	3 (7.40E−03)	
MP:0005344	Increased circulating bilirubin level	4 (7.85E−04)	3 (4.60E−03)	3 (1.69E−02)	
MP:0000598	Abnormal liver morphology	6 (2.15E−02)	4 (4.65E−02)	4 (1.67E−01)	
MP:0002221	Abnormal lymph organ size	5 (7.50E−02)	4 (4.82E−02)	4 (1.84E−01)	
MP:0004952	Increased spleen weight	3 (2.15E−02)	2 (4.82E−02)	2 (1.12E−01)	
MP:0012106	Impaired exercise endurance	2 (3.71E−02)	1 (1.11E−01)	1 (1.92E−01)	
MP:0005584	Abnormal enzyme/coenzyme activity	4 (4.21E−03)	3 (1.11E−02)	3 (3.37E−02)	
MP:0002942	Decreased circulating alanine transaminase level	2 (2.69E−02)	0 (NA)	0 (NA)	
MP:0002943	Abnormal lactate dehydrogenase level	2 (2.32E−02)	1 (9.55E−02)	1 (1.66E−01)	
Phenotypes with more genes in the glucose to acetyl-CoA pathway (aerobic) than glucose to pyruvate	
MP:0030970	Ketosis	1 (1.44E−01)	1 (8.22E−02)	3 (2.23E−03)	
MP:0002575	Increased circulating ketone body level	0 (NA)	0 (NA)	2 (1.74E−02)	
MP:0003059	Decreased insulin secretion	2 (5.96E−02)	2 (3.76E−02)	3 (1.31E−02)	
MP:0005293	Impaired glucose tolerance	3 (6.56E−02)	3 (2.78E−02)	5 (7.40E−03)	
MP:0013403	Abnormal circulating lactate level	1 (1.37E−01)	1 (7.91E−02)	3 (2.23E−03)	
MP:0013405	Increased circulating lactate level	0 (NA)	0 (NA)	2 (1.21E−02)	
MP:0001325	Abnormal retina morphology	2 (3.08E−01)	2 (1.62E−01)	5 (3.20E−02)	
MP:0003731	Abnormal retina outer nuclear layer morphology	2 (1.07E−01)	2 (4.93E−02)	3 (2.50E−02)	
MP:0008931	Abnormal paraventricular thalamic nucleus morphology	0 (NA)	0 (NA)	1 (1.98E−02)	
MP:0004079	Abnormal putamen morphology	0 (NA)	0 (NA)	1 (4.09E−02)	
MP:0020554	Decreased stria medullaris size	0 (NA)	0 (NA)	1 (2.50E−02)	
MP:0009976	Abnormal cerebellar peduncle morphology	0 (NA)	0 (NA)	1 (4.60E−02)	
MP:0003984	Embryonic growth retardation	1 (4.79E−01)	1 (3.24E−01)	4 (3.20E−02)	
Phenotypes with more genes in both the lactate/acetyl-CoA pathways than in glucose to pyruvate	
MP:0000187	Abnormal triglyceride level	3 (1.00E−01)	2 (1.10E−01)	4 (4.09E−02)	
MP:0005318	Decreased triglyceride level	2 (1.38E−01)	1 (2.20E−01)	3 (4.14E−02)	
MP:0002078	Abnormal glucose homeostasis	4 (1.59E−01)	3 (1.16E−01)	7 (2.02E−02)	
MP:0005291	Abnormal glucose tolerance	4 (3.98E−02)	3 (4.65E−02)	6 (7.40E−03)	
MP:0001560	Abnormal circulating insulin level	3 (7.85E−02)	2 (9.49E−02)	4 (2.50E−02)	
MP:0002727	Decreased circulating insulin level	3 (2.69E−02)	2 (5.22E−02)	3 (3.20E−02)	
MP:0002080	Prenatal lethality	5 (1.59E−01)	3 (1.99E−01)	8 (2.60E−02)	
MP:0005266	Abnormal metabolism	6 (4.14E−03)	4 (1.74E−02)	5 (2.50E−02)	
MP:0010768	Mortality/aging	9 (1.57E−01)	7 (7.09E−02)	14 (1.31E−02)	
MP:0010769	Abnormal survival	8 (1.74E−01)	6 (1.00E−01)	13 (1.63E−02)	
MP:0008762	Embryonic lethality	4 (1.74E−01)	3 (1.34E−01)	8 (1.31E−02)	
MP:0010770	Preweaning lethality	8 (1.35E−01)	3 (1.11E−01)	13 (7.40E−03)	
MP:0011100	Preweaning lethality, complete penetrance	4 (1.57E−01)	3 (1.11E−01)	6 (4.60E−02)	
If a parent MP term and its child were annotated to the same numbers of terms, only the child term is shown. For example, “abnormal insulin secretion” (parent—not shown) and “decreased insulin secretion” (child—shown) had 2,2,3 terms.

Phenotypes associated with variant forms of genes of anaerobic metabolism represented abnormalities in the erythropoietic system such as “reticulocytosis” (MP:0002640), “hemolytic anemia” (MP:0001585), “decreased hematocrit” (MP:0000208), and “decreased erythrocyte cell number” (MP:0002875). Phenotypes indirectly associated with the erythropoietic system included “increased circulating bilirubin level” (MP:005344) and “abnormal lactate dehydrogenase level” (MP:002943), 2 clinical features of hemolytic anemia (Dhaliwal et al. 2004). We also observed an increase in the number of genes associated with “impaired exercise endurance” (MP:0012106). One phenotype, “decreased circulating alanine transaminase level” (MP:0002942), was not associated with any genes in the glycolysis pathway, but was associated with 2 genes, Ldha (MGI:96759) and Slc16a1 (MGI:106013) in the downstream conversion of pyruvate to lactate.

Phenotypes associated with variant forms of genes of aerobic metabolism represented abnormalities such as “impaired glucose tolerance” (MP:0005293), “decreased insulin secretion” (MP:0003059), “increased circulating lactate level” (MP:0013405), and “ketosis” (MP:0030970) (a phenotype associated with rapid lipid breakdown). Several phenotypes with central nervous system defects, e.g. “abnormal putamen morphology” (MP:0004079) and “decreased stria medullaris size” (MP:0020544), were associated with genes in the pyruvate to acetyl-CoA portion of the pathway and not the glycolytic portion. Another phenotype, “increased circulating lactate level” (MP:0013405), was only associated with genes in the downstream pyruvate to acetyl-CoA pathway and may indicate a metabolic switch to the pyruvate pathway when the aerobic pathway is not functional.

We also saw general phenotypes associated with genes that function in both the lactate and acetyl-CoA pathways when compared with core glycolysis, e.g. “abnormal triglyceride level” (MP:0000187), abnormal glucose tolerance (MP:0005291), and “decreased circulating insulin level” (MP:0002727).

Phenotypes associated with specific pathways can be used to identify transcription factors causally related to the pathway

Since the GO-CAM models we created represent genes that are causally connected in a pathway, we hypothesized that we could make predictions about other genes that might be causally connected to the same pathways based on their phenotypes. Specifically, we attempted to identify DNA-binding transcription factors causally connected to gluconeogenesis by searching for transcription factors that shared enriched phenotypes. We used the MGI “Genes and Markers” query tool to search for genes annotated with the Gene Ontology term “DNA binding transcription factor activity” (GO:0003700) and the most specific enriched gluconeogenesis phenotypes with FDR < 0.05. The search yielded 93 genes annotated as transcription factors (Table S5), 11 of which were associated with 3 or more phenotypes (Table 4). We then searched for published evidence of causal connections to gluconeogenesis for each of these 11 genes as follows:

Cebpa (MGI:99480) 6 phenotypes: Cebpa directly regulates gluconeogenesis by regulating the gluconeogenic enzymes phosphoenolpyruvate carboxykinase and glucose-6-phosphatase (Pedersen et al. 2007). Furthermore, knocking out Cebpa (Cebpatm1Gjd, MGI:2177053) results in “abnormal gluconeogenesis” (MP:0003383) (Wang et al. 1995).

Bmal1 (MGI:1096381) 5 phenotypes: Bmal1 regulates gluconeogenic gene expression by interaction with Hdac5 (Li et al. 2018). Bmal1 interacts with Clock to form a transcription factor complex (Huang et al. 2012) involved in the circadian regulation of glucose homeostasis (Schiaffino et al. 2016).

Cebpb (MGI:88373) 4 phenotypes: Cebpb directly regulates the expression of the phosphoenolpyruvate carboxykinase gene (Arizmendi et al. 1999).

Pparg (MGI:97747) 4 phenotypes: The Pparg gene is generally thought to have a negative effect on gluconeogenesis and is involved in adipocyte differentiation (Way et al. 2001; Pendse et al. 2010; Hernandez-Quiles et al. 2021). However, Pparg may indirectly upregulate gluconeogenesis through Pomc signaling in response to a high fat diet (Long et al. 2014). The Pparg2 isoform also upregulates pyruvate carboxylase gene expression (Jitrapakdee et al. 2005). Three of the 4 phenotypes associated with Pparg involve triglyceride metabolism. These could all reflect its overall involvement in controlling triglyceride metabolism (Semple et al. 2006), consistent with the fourth “hepatic steatosis” phenotype.

Ncoa3 (MGI:1276535) 3 phenotypes: We were not able to find direct evidence of a role for this gene in gluconeogenesis. It is well documented in the control of lipid metabolism (Astapova 2016), and the members of the Ncoa family are thought to be overall regulators of metabolic processes (York and O’Malley 2010).

Foxa2 (MGI:1347476) and Foxo1 (MGI:1890077) 3 phenotypes each: Foxa2 in concert with Foxo1 upregulates gluconeogenesis and in particular upregulates phosphoenolpyruvate carboxykinase in the fasting liver (Zhang et al. 2005; Puigserver and Rodgers 2006).

Klf15 (MGI:1929988) 3 phenotypes: Klf15 is one of the genes responsible for regulating the switch downregulating lipogenesis and upregulating gluconeogenesis during fasting, upregulating the expression of Pck1 (Takeuchi et al. 2016). Additional results support the role of Klf15 in regulating gluconeogenesis through mutational studies and transcription assays (Teshigawara et al. 2005; Gray et al. 2007).

Nr1h4 (MGI:1352464) 3 phenotypes: Metabolites of sesame oil antagonize Nr1h4 function, resulting in a decrease in mRNA levels encoding phosphoenolpyruvate carboxykinase and glucose 6-phosphatase (Sasaki et al. 2022). Activation of Nr1h4 by cholic acids, however, resulted in a decrease in the expression of gluconeogenic genes, likely by indirectly regulating Foxo1 and Hnf4a (Yamagata et al. 2004). These results may be a reflection of the interaction of Nr1h4 with different partners, resulting in different outcomes (Yamagata et al. 2004; Xu et al. 2018).

Ppara (MGI:104740) 3 phenotypes: Ppara plays a role in regulating many metabolic genes including those that function in gluconeogenesis and lipid metabolism (Preidis et al. 2017). Like other transcription factors we have identified in our query, the exact nature of Ppara's activity may depend on its partners and other transcription factors that are expressed in different physiological states (Kersten 2014).

Clock (MGI:99698) 3 phenotypes: Clock and Bmal1 (above) interact to synchronize gluconeogenesis with circadian rhythms (Schiaffino et al. 2016; Li et al. 2018).

Table 4. Transcription factors associated with 3 or more gluconeogenesis phenotypes.

Number of phenotypes per gene	“Hypoglycemia” (MP:0000189)	“Increased liver triglyceride level” (MP:0009355)	“Decreased circulating triglyceride level” (MP:0002644)	“Increased liver glycogen level” (MP:0010400)	“Decreased granulocyte number” (MP:0000334)	“Abnormal gluconeogenesis” (MP:0003383)	“Increased circulating triglyceride level” (MP:0001552)	“Hepatic steatosis” (MP:0002628)	
6	Cebpa	Cebpa			Cebpa	Cebpa	Cebpa	Cebpa	
5		Bmal1	Bmal1		Bmal1	Bmal1	Bmal1		
4	Cebpb		Cebpb			Cebpb	Cebpb		
4		Pparg	Pparg				Pparg	Pparg	
3		Ncoa3	Ncoa3	Ncoa3					
3	Foxa2	Foxa2					Foxa2		
3			Foxo1			Foxo1	Foxo1		
3	Klf15					Klf15	Klf15		
3		Nr1h4					Nr1h4	Nr1h4	
3	Ppara						Ppara	Ppara	
3						Clock	Clock	Clock	

Although it is standard practice at MGI to infer the biological processes in which a gene is involved using the “inferred from mutant phenotype” evidence code, GO and phenotype annotation are done separately. Of the 11 genes that we identified through our phenotypic analysis, only 2, Ppara and Foxo1, were already annotated in the separate GO curation process as being involved in gluconeogenesis or its regulation. Six of the 11 transcription factors we identified (Cebpa, Bmal1, Cebpb, Foxo1, Klf15, Clock) were also annotated to “abnormal gluconeogenesis” (MP:0003383). To test whether the transcription factors were preferentially annotated to gluconeogenic phenotypes, we ran a VLAD analysis on the transcription factors alone and compared with the gluconeogenic genes. As expected, the top scoring phenotypes were related to the pre-chosen gluconeogenic phenotypes. However, the transcription factors were also significantly enriched for other phenotypes that were not enriched in the gluconeogenic gene analysis, e.g. “abnormal ovulation cycle” (MP:0009344) and “abnormal hair cycle anagen phase” (MP:0008858) (Supplementary Table S6, analysis date 2023 July 3). These results imply that the transcription factors we identified were not a priori or specifically limited to the study of gluconeogenesis.

Discussion

Here, we have shown that human GO-CAM models of glycolysis, gluconeogenesis, and pyruvate metabolism serve as accurate templates for modeling the orthologous mouse pathways, that VLAD analysis of phenotypes associated with mutant mouse genes associated distinct but related phenotypes with mutational disruptions of genes involved in each pathway, and that analysis of common enriched phenotypes further supported the identification of transcription factors that play major roles in regulating these pathways. These results are not surprising: These pathways, their transcriptional regulation, and their pathophysiology have been very extensively characterized experimentally in both species and shown to be closely related. Neither GO-CAM construction nor VLAD analysis was trained or optimized for these domains of biology, however, so these results are useful as a test of a general strategy for pathway annotation and analysis, and suggest that it can be extended to less exhaustively studied processes and less well characterized species.

Reactome GO-CAMs can bridge the biology of model organisms and humans

Reactome annotations of human pathways of glycolysis, gluconeogenesis, and pyruvate metabolism to lactate or acetyl-CoA identified functions for 60 human gene products, 44 of which have high-confidence mouse structural orthologs (Table 1), enabling the streamlined generation of mouse GO-CAM models by swapping human gene products for their mouse counterparts. This successful manual pathway construction exercise is a proof of principle for automating the construction of model organism template pathways useful as starting points for expert manual curation and addition of organism-specific supporting experimental evidence. Indeed, Reactome already computationally predicts pathways for 15 model organisms using orthologs of human proteins to systematically populate reactions for these other species, a process that could be extended to any model organism whose proteome is known (https://reactome.org/documentation/inferred-events).

Phenotype enrichment analysis of pathway-specific gene lists identifies and distinguishes pathway-specific phenotypes

We separated the genes annotated in our mouse GO-CAM models into 3 groups as follows: specific for glycolysis, specific for gluconeogenesis, and shared genes. Phenotype enrichment on the 3 sets showed that the enriched phenotypes fell into different classes for each set (Figs. 2 and 3). Phenotypes associated with alleles of the genes shared by glycolysis and gluconeogenesis compared with genes specific to each pathway were overrepresented in an early embryonic lethal phenotype, reflecting the central role of glucose metabolism in mammalian physiology.

Variants of the glycolysis-specific genes, compared to gluconeogenesis-specific or shared genes, were more significantly enriched for phenotypes associated with abnormalities of red blood cells, consistent with the absolute dependence of these cells on anaerobic glycolysis for ATP.

Variants of gluconeogenesis-specific genes were overrepresented in the abnormal phenotypes of hypoglycemia, increased triglyceride levels and abnormal kidney and liver glycogen levels. Increased glycogen levels in the kidney and liver are explained by the identities of the genetic defects associated with the phenotypes: mutations that inactivate Slc37a4 (MGI:1316650) and G6pc (MGI:95607). The products of these 2 genes normally mediate the transport and dephosphorylation of glucose-6-phosphate newly synthesized by gluconeogenesis, enabling release of free glucose into the circulation–inactivation of either leads directly to hypoglycemia (Froissart et al. 2011). Abnormal accumulation of intracellular glucose-6-phosphate drives glycogen synthesis by mass action (Adeva-Andany et al. 2016). “Increased liver triglyceride level” was associated with mutations of 3 genes in the gluconeogenesis-specific gene set, G6pc (MGI:95607), Pck1 (MGI:97501), and Slc25a13 (MGI:1354721), and is consistent with phenotypes associated with excess accumulation of glucose-6-phosphate observed in a G6pc knockout mouse and in human patients with partial loss of G6PC function (Hoogerland et al. 2021).

To test the hypothesis that phenotypes could be further distinguished by joining analysis of glycolysis with analysis of the aerobic and anaerobic fates of its major product, pyruvate, we performed phenotype enrichment analyses on 3 sets of genes as follows: those specific for canonical glycolysis (conversion of glucose to pyruvate), those for canonical glycolysis plus downstream anaerobic metabolism of pyruvate to lactate, and those specific for canonical glycolysis plus downstream aerobic metabolism of pyruvate to acetyl-CoA (Table 3). Analysis of the glycolysis + lactate/anaerobic gene set revealed additional associations of genes with red blood cell phenotypes. These phenotypes reflect the dependence of red blood cells on anaerobic glycolysis to lactate as an energy source (Goto et al. 2019); (Rose and Warms 1966); (TeSlaa et al. 2021). A non-erythropoietic phenotype identified in this analysis was “impaired exercise endurance”, consistent with enhanced anaerobic glucose metabolism during intense exercise (Alberti 1977; Robergs et al. 2004). When we add the genes for the acetyl-CoA/aerobic pathway, additional associations with glucose homeostasis phenotypes are detected. These associations are in concordance with known physiology. Glycolysis and the downstream aerobic metabolism of glucose through pyruvate conversion to acetyl-CoA and the tricarboxylic acid cycle are key steps in the regulation of insulin secretion in response to glucose in pancreatic beta cells (Newgard and McGarry 1995).These results show that we can identify connected pathways that give rise to distinct phenotypes and to distinguish subnetworks in a branching system that have discreet biological outcomes.

Enriched phenotypes can be used to identify genes that are causally connected to a pathway

To test our ability to identify genes that are not conventionally associated with a pathway but that are causally connected to it by shared phenotypes, we used enriched phenotypes in the gluconeogenesis pathway to search for transcription factor genes whose variants were associated with these same phenotypes. Bmal1 (Li et al. 2018) Cebpa (Pedersen et al. 2007), Cebpb (Arizmendi et al. 1999), and Foxa2 (Zhang et al. 2005; Puigserver and Rodgers 2006), either have a direct or indirect effect on the expression of gluconeogenesis genes. Several genes, such as Ncoa3 (York and O’Malley 2010) and Pparg (Semple et al. 2006), are also connected in less direct ways to gluconeogenesis and may be connected to other pathways causally connected to gluconeogenesis such as glucagon signaling and lipid metabolism. Some genes, like Ppara, regulate gluconeogenesis as well as other related, connected metabolic processes like lipid metabolism (Rakhshandehroo et al. 2007). Our results are consistent with a complex regulation of gluconeogenesis where a variety of transcription factors act together both directly and indirectly to coordinate metabolism depending on specific contexts (Pei et al. 2006; Li et al. 2018; Qu et al. 2021; De Sousa-Coelho et al. 2023).

In an extension of previous work demonstrating that phenotypes can be used to predict GO annotations, our work further shows the integrated nature of phenotypic analysis and functional inference (Ascensao et al. 2014). With the exception of 2 transcription factors, Ppara and Foxo1, none of the genes we identified by shared phenotypes were yet annotated to gluconeogenesis GO terms even though 6 of them were annotated to the phenotype term “abnormal gluconeogenesis” (MP:0003383). Although it is serendipitous that some of the transcription factors had not yet been annotated to GO based on phenotypes associated with gluconeogenesis, this suggests that genes were not “selected” with a preexisting annotation focus on gluconeogenesis. Instead, it shows that phenotypes can be used to inform us about gene function and that gene function may be useful for inferring phenotypes.

We hypothesize that as we integrate models of additional regulatory processes with models of metabolism, we will be able to discriminate among the networks that these transcription factors modulate and identify specific gene targets within them. Concatenation of regulatory and housekeeping processes will enable construction of vast network hairballs; the limited study described here suggests that it may be possible to manage and interpret them to generate integrated models of cell physiology useful for modeling and experimental design.

Limitations of our current approach and future directions

To increase the sensitivity of our search for abnormal phenotypes associated with pathways, we have grouped the effects of all paralogs that normally mediate each activity in glycolysis, gluconeogenesis, and pyruvate metabolism. Despite this grouping, known tissue- and process-specific effects, such as the impact of glycolysis defects on red cell function and sperm cell maturation and function, and of gluconeogenesis defects on glucose homeostasis can be detected above background significance levels (Table S1). Nevertheless, we have missed other known effects, such as the opposing effects of variant forms of GCK (human)/Gck (mouse—MGI:1270854) that have abnormally high and low affinities for glucose on insulin secretion in response to varying blood glucose levels. This result suggests that our analysis could be refined further by taking into account the nature of the alleles that result in mutant phenotypes, the background strains on which alleles are present, and the physiological conditions in which phenotypes are measured. We also failed to observe red blood cell phenotypes for all genes required for conversion of glucose to pyruvate. These failures have several possible explanations. First, the collection of known genetic variants is far from saturating, and phenotypic characterization of existing variants has not been exhaustive. Second, many variants in these genes are likely to be incompatible with early embryonic development, so their effects on postnatal functions would not be detectable. Third, most of the steps of glycolysis are mediated by sets of paralogous genes that differ from one another in their regulation and their tissue-specific expression, but it is plausible that phenotypes due to variants in 1 member of a set could be masked by even low level expression of other (wild-type) set members.

As in all bioinformatics approaches, our analyses are limited by the data available. We chose well-studied pathways for this work to test whether when the data exist, we can use those data to generate specific conclusions. In cases where phenotypic data are not available or comprehensive for a large number of genes in a pathway or the pathway is not well-defined, our approach will be limited due to lack of enrichment or multiple genes showing the same phenotype and phenotypes may be missed. There will be a “tipping point” in which enough genes in a pathway will be annotated to a phenotype to result in a significant enrichment or will give a clear indication that disturbing the pathway results in a phenotype. However, once the “tipping point” is reached, we can predict that perturbing the other genes in the causal chain will result in the same phenotype.

Relating biochemical pathways to phenotypes associated with specific genes provides a tool to predict and test the effects of genetic or pharmacological manipulation on the physiology of an organism. Identifying the elements of pathways, including ones with complicated structures, that result in a phenotype should better support targeted approaches to modify the pathways and change phenotypic outcomes. Beyond the classic approach of providing exogenous sources of key molecules missing due to a gene defect, we can envision strategies to manipulate pathways that are hypermorphic and are deleterious to an organism, with pharmacologic inhibitors or RNAi methods. Finally, these types of analyses may help us to dissect the complex interactions of drugs with biochemical pathways to better understand mechanisms and potential side effects that could be addressed with co-therapies. For example, metformin is one of the primary pharmaceuticals used to treat type 2 diabetes. Metformin lowers gluconeogenesis in the liver, but its actions are complex and not completely understood (Rena et al. 2017; Flory and Lipska 2019). As we further enrich our representation of glucose metabolism to include metabolic and regulatory networks, it might be possible to perform a detailed analysis of the effects (phenotypes) of metformin treatment in mice and compare that to network analysis similar to the ones we performed in this work. The comparison can be used to identify potential branches in the network where metformin acts by looking for shared or unique phenotypes.

The available mouse phenotyping data themselves present some limitations. Different mutations in the same gene can have variable phenotypes due to the molecular nature of the mutation (e.g. null or knockout vs point mutation, or isoform-specific mutations). Mutations may be hypomorphic, gain or loss of function or even silent in the absence of an inducer or stressor. Phenotypic variability due to strain-specific variants in the mouse experimental genetic background can contribute to conflicting results when the data is inferred at the gene level (Perry et al. 2020). Genetic sex, microbiota, and parental origin (imprinting) can also influence phenotypes. Outliers in the data should be assessed for these potential causes in variability.

In summary, we have shown that we can convert a set of human biochemical pathways to corresponding pathways in the laboratory mouse. We can then take advantage of the rich genetic work that has been done in the mouse to analyze and determine the phenotypic consequences of perturbations of genes in those pathways. In this way, we integrate the information about human biology from Reactome with model organism biology from MGI. We envision that this approach can be extended to other metabolic processes and other model organisms and the study of pathway-phenotype connections can be used not only to understand the similarities of the pathways but as a testing ground for manipulation of pathways in more experimentally tractable organisms than human.

Supplementary Material

iyad152_Supplementary_Data

Acknowledgments

We dedicate this paper to the memory of Michael Ashburner. The work described here comes out of projects he helped to initiate and guide, and we benefitted immeasurably from his support.

Data availability

Full results of our analyses are given in Tables 1–4 and in the Supplementary Files for this paper. The full Reactome and MGI data sets are freely available to all users at www.reactome.org and https://www.informatics.jax.org/, respectively.

Supplemental material available at GENETICS online.

Funding

This work was supported by NIH grants U41 HG002273 and U24 HG012212 (GO Consortium), U24 HG012198 (Reactome), U24 HG002223 (WormBase), U41 HG000330 (The Mouse Genome Database), and U24 HG011851 (Pathways2GO).
==== Refs
Literature cited

Adeva-Andany  MM, González-Lucán  M, Donapetry-García  C, Fernández-Fernández  C, Ameneiros-Rodríguez  E. Glycogen metabolism in humans. BBA Clin. 2016;5 :85–100. doi:10.1016/j.bbacli.2016.02.001.27051594
Alberti  KG . The biochemical consequences of hypoxia. J Clin Pathol Suppl (R Coll Pathol). 1977;11 (1 ):14–20. doi:10.1136/jcp.s3-11.1.14.198434
Arizmendi  C, Liu  S, Croniger  C, Poli  V, Friedman  JE. The transcription factor CCAAT/enhancer-binding protein beta regulates gluconeogenesis and phosphoenolpyruvate carboxykinase (GTP) gene transcription during diabetes. J Biol Chem. 1999;274 (19 ):13033–13040. doi:10.1074/jbc.274.19.13033.10224054
Ascensao  JA, Dolan  ME, Hill  DP, Blake  JA. Methodology for the inference of gene function from phenotype data. BMC Bioinformatics. 2014;15 (1 ):405. doi:10.1186/s12859-014-0405-z.25495798
Astapova  I . Role of co-regulators in metabolic and transcriptional actions of thyroid hormone. J Mol Endocrinol. 2016;56 (3 ):73–97. doi:10.1530/JME-15-0246.26673411
Balakrishnan  R, Harris  MA, Huntley  R, Van Auken  K, Cherry  JM. A guide to best practices for Gene Ontology (GO) manual annotation. Database (Oxford). 2013;2013 :bat054. doi:10.1093/database/bat054.23842463
Barton  AR, Hujoel  MLA, Mukamel  RE, Sherman  MA, Loh  P-R. A spectrum of recessiveness among Mendelian disease variants in UK Biobank. Am J Hum Genet. 2022;109 (7 ):1298–1307. doi:10.1016/j.ajhg.2022.05.008.35649421
Bello  SM, Eppig  JT, the MGI Software Group. Inferring gene-to-phenotype and gene-to-disease relationships at Mouse Genome Informatics: challenges and solutions. J Biomed Semant. 2016;7 (1 ):14. doi:10.1186/s13326-016-0054-4.
Bello  SM, Perry  MN, Smith  CL. Know your model: a brief history of making mutant mouse genetic models. Lab Anim (NY). 2021;50 (10 ):263–266. doi:10.1038/s41684-021-00853-5.34561680
Bello  SM, Smith  CL, Eppig  JT. Allele, phenotype and disease data at Mouse Genome Informatics: improving access and analysis. Mamm Genome. 2015;26 (7–8 ):285–294. doi:10.1007/s00335-015-9582-y.26162703
Blake  JA, Baldarelli  R, Kadin  JA, Richardson  JE, Smith  CL, Bult  CJ, Anagnostopoulos  AV, Beal  JS, Bello  SM, Blodgett  O, et al  Mouse Genome Database (MGD): knowledgebase for mouse–human comparative biology. Nucleic Acids Res. 2021;49 (D1 ):D981–D987. doi:10.1093/nar/gkaa1083.33231642
Bulusu  V, Prior  N, Snaebjornsson  MT, Kuehne  A, Sonnen  KF, Kress  J, Stein  F, Schultz  C, Sauer  U, Aulehla  A. Spatiotemporal analysis of a glycolytic activity gradient linked to mouse embryo mesoderm development. Dev Cell. 2017;40 (4 ):331–341.e4. doi:10.1016/j.devcel.2017.01.015.28245920
Chourpiliadis  C, Mohiuddin  SS. Biochemistry, gluconeogenesis. In: StatPearls. Treasure Island (FL): StatPearls Publishing; 2022.
Cori  CF . The glucose-lactic acid cycle and gluconeogenesis. Curr Top Cell Regul. 1981;18 :377–387. doi:10.1016/B978-0-12-152818-8.50028-1.7273846
De Sousa-Coelho  AL, Gacias  M, O’Neill  BT, Relat  J, Link  W, Haro  D, Marrero  PF. FOXO1 represses PPARα-mediated induction of FGF21 gene expression. Biochem Biophys Res Commun. 2023;644 :122–129. doi:10.1016/j.bbrc.2023.01.012.36640666
de Visser  JAGM, Cooper  TF, Elena  SF. The causes of epistasis. Proc Biol Sci. 2011;278 (1725 ):3617–3624. doi:10.1098/rspb.2011.1537.21976687
Dhaliwal  G, Cornett  PA, Tierney  LM. Hemolytic anemia. Am Fam Physician. 2004;69 :2599–2606.15202694
Felig  P . The glucose-alanine cycle. Metabolism. 1973;22 (2 ):179–207. doi:10.1016/0026-0495(73)90269-2.4567003
Flory  J, Lipska  K. Metformin in 2019. JAMA. 2019;321 (19 ):1926–1927. doi:10.1001/jama.2019.3805.31009043
Froissart  R, Piraud  M, Boudjemline  AM, Vianey-Saban  C, Petit  F, Hubert-Buron  A, Eberschweiler  PT, Gajdos  V, Labrune  P. Glucose-6-phosphatase deficiency. Orphanet J Rare Dis. 2011;6 (1 ):27. doi:10.1186/1750-1172-6-27.21599942
Gillespie  M, Jassal  B, Stephan  R, Milacic  M, Rothfels  K, Senff-Ribeiro  A, Griss  J, Sevilla  C, Matthews  L, Gong  C, et al  The reactome pathway knowledgebase 2022. Nucleic Acids Res. 2022;50 (D1 ):D687–D692. doi:10.1093/nar/gkab1028.34788843
Good  BM, Van Auken  K, Hill  DP, Mi  H, Carbon  S, Balhoff  JP, Albou  L-Philippe, Thomas  PD, Mungall  CJ, Blake  JA, et al  Reactome and the Gene Ontology: digital convergence of data resources. Bioinformatics. 2021;37 (19 ):3343–3348. doi:10.1093/bioinformatics/btab325.33964129
Goto  T, Ubukawa  K, Kobayashi  I, Sugawara  K, Asanuma  K, Sasaki  Y, Guo  Y-M, Takahashi  N, Sawada  K, Wakui  H, et al  ATP produced by anaerobic glycolysis is essential for enucleation of human erythroblasts. Exp Hematol. 2019;72 :14–26.e1. doi:10.1016/j.exphem.2019.02.004.30797950
Gray  S, Wang  B, Orihuela  Y, Hong  E-G, Fisch  S, Haldar  S, Cline  GW, Kim  JK, Peroni  OD, Kahn  BB, et al  Regulation of gluconeogenesis by Krüppel-like factor 15. Cell Metab. 2007;5 (4 ):305–312. doi:10.1016/j.cmet.2007.03.002.17403374
Groza  T, Gomez  FL, Mashhadi  HH, Muñoz-Fuentes  V, Gunes  O, Gunes  O, Wilson  R, Cacheiro  P, Frost  A, Keskivali-Bond  P, et al  The International Mouse Phenotyping Consortium: comprehensive knockout phenotyping underpinning the study of human disease. Nucleic Acids Res. 2023;51 (D1 ):D1038–D1045. doi:10.1093/nar/gkac972.36305825
Hernandez-Quiles  M, Broekema  MF, Kalkhoven  E. PPARgamma in metabolism. Immunity, and cancer: unified and diverse mechanisms of action. Front Endocrinol (Lausanne). 2021;12 :624112. doi:10.3389/fendo.2021.624112.33716977
Hill  DP, D’Eustachio  P, Berardini  TZ, Mungall  CJ, Renedo  N, Blake  JA. Modeling biochemical pathways in the gene ontology. Database. 2016;2016 :baw126. doi:10.1093/database/baw126.27589964
Hoogerland  JA, Peeks  F, Hijmans  BS, Wolters  JS, Kooijman  S, Bos  T, Bleeker  A, Dijk  TH, Wolters  H, Gerding  A, et al  Impaired very-low-density lipoprotein catabolism links hypoglycemia to hypertriglyceridemia in glycogen storage disease type Ia. J Inherit Metab Dis. 2021;44 (4 ):879–892. doi:10.1002/jimd.12380.33739445
Huang  N, Chelliah  Y, Shan  Y, Taylor  CA, Yoo  S-H, Partch  C, Green  CB, Zhang  H, Takahashi  JS. Crystal structure of the heterodimeric CLOCK:BMAL1 transcriptional activator complex. Science. 2012;337 (6091 ):189–194. doi:10.1126/science.1222804.22653727
Jitrapakdee  S, Slawik  M, Medina-Gomez  G, Campbell  M, Wallace  JC, Sethi  JK, O'Rahilly  S, Vidal-Puig  AJ. The peroxisome proliferator-activated receptor-gamma regulates murine pyruvate carboxylase gene expression in vivo and in vitro. J Biol Chem. 2005;280 (29 ):27466–27476. doi:10.1074/jbc.M503836200.15917242
Kersten  S . Integrated physiology and systems biology of PPARα. Mol Metab. 2014;3 (4 ):354–371. doi:10.1016/j.molmet.2014.02.002.24944896
Kingdom  R, Tuke  M, Wood  A, Beaumont  RN, Frayling  TM, Weedon  MN, Wright  CF. Rare genetic variants in genes and loci linked to dominant monogenic developmental disorders cause milder related phenotypes in the general population. Am J Hum Genet. 2022;109 (7 ):1308–1316. doi:10.1016/j.ajhg.2022.05.011.35700724
Kingdom  R, Wright  CF. Incomplete penetrance and variable expressivity: from clinical studies to population cohorts. Front Genet. 2022;13 :920390. doi:10.3389/fgene.2022.920390.35983412
Le Novère  N, Hucka  M, Mi  H, Moodie  S, Schreiber  F, Sorokin  A, Demir  E, Wegner  K, Aladjem  MI, Wimalaratne  SM, et al  The systems biology graphical notation. Nat Biotechnol. 2009;27 (8 ):735–741. doi:10.1038/nbt.1558.19668183
Li  J, Lv  S, Qiu  X, Yu  J, Jiang  J, Jin  Y, Guo  W, Zhao  R, Zhang  Z-N, Zhang  C, et al  BMAL1 functions as a cAMP-responsive coactivator of HDAC5 to regulate hepatic gluconeogenesis. Protein Cell. 2018;9 (11 ):976–980. doi:10.1007/s13238-018-0514-y.29508277
Long  L, Toda  C, Jeong  JK, Horvath  TL, Diano  S. PPARγ ablation sensitizes proopiomelanocortin neurons to leptin during high-fat feeding. J Clin Invest. 2014;124 (9 ):4017–4027. doi:10.1172/JCI76220.25083994
Mirshahi  UL, Colclough  K, Wright  CF, Wood  AR, Beaumont  RN, Tyrrell  J, Laver  TW, Stahl  R, Golden  A, Goehringer  JM, et al  Reduced penetrance of MODY-associated HNF1A/HNF4A variants but not GCK variants in clinically unselected cohorts. Am J Hum Genet. 2022;109 (11 ):2018–2028. doi:10.1016/j.ajhg.2022.09.014.36257325
Newgard  CB, McGarry  JD. Metabolic coupling factors in pancreatic beta-cell signal transduction. Annu Rev Biochem. 1995;64 (1 ):689–719. doi:10.1146/annurev.bi.64.070195.003353.7574498
Pedersen  TA, Bereshchenko  O, Garcia-Silva  S, Ermakova  O, Kurz  E, Mandrup  S, Porse  BT, Nerlov  Cs. Distinct C/EBPalpha motifs regulate lipogenic and gluconeogenic gene expression in vivo. EMBO J. 2007;26 (4 ):1081–1093. doi:10.1038/sj.emboj.7601563.17290224
Pei  L, Waki  H, Vaitheesvaran  B, Wilpitz  DC, Kurland  IJ, Tontonoz  P. NR4A orphan nuclear receptors are transcriptional regulators of hepatic glucose metabolism. Nat Med. 2006;12 (9 ):1048–1055. doi:10.1038/nm1471.16906154
Pendse  AA, Johnson  LA, Tsai  YS, Maeda  N. Pparg-P465L mutation worsens hyperglycemia in Ins2-Akita female mice via adipose-specific insulin resistance and storage dysfunction. Diabetes. 2010;59 (11 ):2890–2897. doi:10.2337/db10-0673.20724579
Perry  MN, Bello  SM, Smith  CL. Know your model: why mouse inbred strain contribution matters. Lab Anim (NY). 2020;49 (5 ):133–134. doi:10.1038/s41684-020-0525-1.32284572
Preidis  GA, Kim  KH, Moore  DD. Nutrient-sensing nuclear receptors PPARα and FXR control liver energy balance. J Clin Invest. 2017;127 (4 ):1193–1201. doi:10.1172/JCI88893.28287408
Prochownik  EV, Wang  H. The metabolic fates of pyruvate in normal and neoplastic cells. Cells. 2021;10 (4 ):762. doi:10.3390/cells10040762.33808495
Puigserver  P, Rodgers  JT. Foxa2, a novel transcriptional regulator of insulin sensitivity. Nat Med. 2006;12 (1 ):38–39. doi:10.1038/nm0106-38.16397559
Qu  M, Qu  H, Jia  Z, Kay  SA. HNF4A defines tissue-specific circadian rhythms by beaconing BMAL1::CLOCK chromatin binding and shaping the rhythmic chromatin landscape. Nat Commun. 2021;12 (1 ):6350. doi:10.1038/s41467-021-26567-3.34732735
Rakhshandehroo  M, Sanderson  LM, Matilainen  M, Stienstra  R, Carlberg  C, de Groot  PJ, Müller  M, Kersten  S. Comprehensive analysis of PPARalpha-dependent regulation of hepatic lipid metabolism by expression profiling. PPAR Res. 2007;2007 :26839. doi:10.1155/2007/26839.18288265
Rena  G, Hardie  DG, Pearson  ER. The mechanisms of action of metformin. Diabetologia. 2017;60 (9 ):1577–1585. doi:10.1007/s00125-017-4342-z.28776086
Richardson  JE, Bult  CJ. Visual annotation display (VLAD): a tool for finding functional themes in lists of genes. Mamm Genome. 2015;26 (9–10 ):567–573. doi:10.1007/s00335-015-9570-2.26047590
Ringwald  M, Richardson  JE, Baldarelli  RM, Blake  JA, Kadin  JA, Smith  C, Bult  CJ. Mouse Genome Informatics (MGI): latest news from MGD and GXD. Mamm Genome. 2022;33 (1 ):4–18. doi:10.1007/s00335-021-09921-0.34698891
Robergs  RA, Ghiasvand  F, Parker  D. Biochemistry of exercise-induced metabolic acidosis. Am J Physiol Regul Integr Comp Physiol. 2004;287 (3 ):R502–R516. doi:10.1152/ajpregu.00114.2004.15308499
Rose  IA, Warms  JV. Control of glycolysis in the human red blood cell. J Biol Chem. 1966;241 (21 ):4848–4854. doi:10.1016/S0021-9258(18)99643-2.4288723
Sasaki  T, Okuda  M, Hong  TW, Watanabe  Y, Takahashi  Y, Shimizu  M, Yamauchi  Y, Sato  R. Sesamin and hepatic metabolites derived from sesamin and episesamin antagonize farnesoid X receptor and reduce the expression of gluconeogenesis-related genes. J Nutr Sci Vitaminol (Tokyo). 2022;68 (1 ):55–64. doi:10.3177/jnsv.68.55.35228496
Schiaffino  S, Blaauw  B, Dyar  KA. The functional significance of the skeletal muscle clock: lessons from Bmal1 knockout models. Skelet Muscle. 2016;6 (1 ):33. doi:10.1186/s13395-016-0107-5.27752300
Semple  RK, Chatterjee  VKK, O’Rahilly  S. PPAR gamma and human metabolic disease. J Clin Invest. 2006;116 (3 ):581–589. doi:10.1172/JCI28003.16511590
Smith  CL, Eppig  JT. The Mammalian Phenotype Ontology as a unifying standard for experimental and high-throughput phenotyping data. Mamm Genome. 2012;23 (9–10 ):653–668. doi:10.1007/s00335-012-9421-3.22961259
Smith  CL, Goldsmith  CA, Eppig  JT. The Mammalian Phenotype Ontology as a tool for annotating, analyzing and comparing phenotypic information. Genome Biol. 2005;6 (1 ):R7. doi:10.1186/gb-2004-6-1-r7.15642099
Takeuchi  Y, Yahagi  N, Aita  Y, Murayama  Y, Sawada  Y, Piao  X, Toya  N, Oya  Y, Shikama  A, Takarada  A, et al  KLF15 enables rapid switching between lipogenesis and gluconeogenesis during fasting. Cell Rep. 2016;16 (9 ):2373–2386. doi:10.1016/j.celrep.2016.07.069.27545894
Teshigawara  K, Ogawa  W, Mori  T, Matsuki  TY, Watanabe  E, Hiramatsu  R, Inoue  H, Miyake  K, Sakaue  H, Kasuga  M. Role of Krüppel-like factor 15 in PEPCK gene expression in the liver. Biochem Biophys Res Commun. 2005;327 (3 ):920–926. doi:10.1016/j.bbrc.2004.12.096.15649433
TeSlaa  T, Bartman  CR, Jankowski  CSR, Zhang  Z, Xu  X, Xu  X, Xing  X, Wang  L, Lu  W, Hui  S, et al  The source of glycolytic intermediates in mammalian tissues. Cell Metab. 2021;33 (2 ):367–378.e5. doi:10.1016/j.cmet.2020.12.020.33472024
The Gene Ontology Consortium . Expansion of the Gene Ontology knowledgebase and resources. Nucleic Acids Res. 2017;45 (D1 ):D331–D338. doi:10.1093/nar/gkw1108.27899567
Thomas  PD, Hill  DP, Mi  H, Osumi-Sutherland  D, Van Auken  K, Carbon  S, Balhoff  JP, Albou  L-P, Good  B, Gaudet  P, et al  Gene Ontology Causal Activity Modeling (GO-CAM) moves beyond GO annotations to structured descriptions of biological functions and systems. Nat Genet. 2019;51 (10 ):1429–1433. doi:10.1038/s41588-019-0500-1.31548717
Wang  ND, Finegold  MJ, Bradley  A, Ou  CN, Abdelsayed  SV, Wilde  Margaret D., Taylor  LR, Wilson  DR, Darlington  GJ. Impaired energy homeostasis in C/EBP alpha knockout mice. Science. 1995;269 (5227 ):1108–1112. doi:10.1126/science.7652557.7652557
Way  JM, Harrington  WW, Brown  KK, Gottschalk  WK, Sundseth  SS, Mansfield  TA, Ramachandran  RK, Willson  TM, Kliewer  SA. Comprehensive messenger ribonucleic acid profiling reveals that peroxisome proliferator-activated receptor gamma activation has coordinate effects on gene expression in multiple insulin-sensitive tissues. Endocrinology. 2001;142 (3 ):1269–1277. doi:10.1210/endo.142.3.8037.11181544
Xu  X, Shi  X, Chen  Y, Zhou  T, Wang  J, Xu  X, Chen  L, Hu  L, Shen  X. HS218 as an FXR antagonist suppresses gluconeogenesis by inhibiting FXR binding to PGC-1α promoter. Metab Clin Exp.  2018;85 :126–138. doi:10.1016/j.metabol.2018.03.016.29577938
Yamagata  K, Daitoku  H, Shimamoto  Y, Matsuzaki  H, Hirota  K, Hirota  K, Ishida  J, Fukamizu  A. Bile acids regulate gluconeogenic gene expression via small heterodimer partner-mediated repression of hepatocyte nuclear factor 4 and Foxo1. J Biol Chem. 2004;279 (22 ):23158–23165. doi:10.1074/jbc.M314322200.15047713
York  B, O’Malley  BW. Steroid receptor coactivator (SRC) family: masters of systems biology. J Biol Chem. 2010;285 (50 ):38743–38750. doi:10.1074/jbc.R110.193367.20956538
Zhang  L, Rubins  NE, Ahima  RS, Greenbaum  LE, Kaestner  KH. Foxa2 integrates the transcriptional response of the hepatocyte to fasting. Cell Metab. 2005;2 (2 ):141–148. doi:10.1016/j.cmet.2005.07.002.16098831
