A framework for application of metabolic modeling in yeast to predict the effects of nsSNV in human orthologs

Reciprocal best BLAST hits method to identify human-yeast ortholog pairs

Identification of candidate enzymes for in vivo validation

We must also consider the case when the product of an enzymatic reaction can be synthesized by alternative means. For this paper, we decided to consider these cases due to literature supporting biochemical importance of dysfunction despite multiple pathways to product.

Additionally, we want to ensure the variable site of interest, the active site, is conserved across ortholog pairs. If functional residues are not conserved, we cannot say that experimental observations resulting from a variation in a yeast enzyme is at all indicative of outcomes of variation in the corresponding human enzyme. This was checked by manual examination of alignments (Figure 1). The minimum requirement for active site conservation is that the residue(s) annotated as the active site residue(s) are identical across species, although many of the human-yeast pairs have blocks of up 10 amino acids conserved across species. 10 of the 11 ortholog pairs were conserved at the active site.

Finally, we want to make sure the ortholog pairs are involved in similar pathways. If a human protein has diverged to act on different metabolites we again cannot use its yeast ortholog as a model to understand effects of variation in human proteins. Pathway information for the human proteins was previously retrieved and summarized in Table 2. Pathway information for the yeast proteins was similarly retrieved from KEGG and recorded. Pathway conservation is indicated by similar pathway map with identical substrates and products. All 6 ortholog pairs were found to participate in the same pathways, acting on the same substrates and yielding the same products. Pathway information for these 6 proteins is included in Table 4. Further support for functional similarity can be seen in Figure 2, showing evidenced interactions (collected and curated following guidance by Lim et. al. [37] with other proteins and the distribution of functional human-yeast orthologs predicted by Isobase [38] among these subsets. A summary schema of the overall method used in identification of these 6 candidate proteins is presented in Figure 3.

Human UniProt Ac	Gene name	Pathways	Yeast UniProt Ac	Yeast gene name	Pathway match?	Details
P13716	ALAD	Porphyrin and clorophyll metabolism	P05373	HEM2	Yes	Same
		Biosynthesis of secondary metabolites			No	Pathway only annotated in yeast; same substrate/product
P22830	FECH	Porphyrin and clorophyll metabolism	P16622	HEM15	Yes	Same
		Biosynthesis of secondary metabolites			No	Pathway only annotated in yeast; same substrate/product
P48637	GSS	Glutathione metabolism	Q08220	GSH2	Yes	Same
P49247	RPIA	Pentose phosphate pathway	Q12189	RKI1	Yes	Same
		Biosynthesis of secondary metabolites			No	Pathway only annotated in yeast; same substrate/product
Q96GX9	APIP	Cysteine and methionine metabolism	P47095	MDE1	Yes	Same
Q9Y2Z4	YARS2	Aminoacyl-tRNA biosynthesis	P48527	MSY1	Yes	Same

This table summarizes the pathway involvement of the six identified candidate proteins. For all 6 proteins, the primary pathway involvement is the same, down to the same products and substrates and the same place within a pathway. For three of the yeast proteins, however, the pathway annotation in KEGG listed an additional pathway not annotated for humans: “Biosynthesis of secondary metabolites.” This is likely due to differential annotation between yeast and humans within KEGG. Because the substrates and products annotated for these pathways were identical to those of the primary pathways, these three proteins were not excluded and went onto analysis in the yeast FBA modeling.

Yeast 7.0 COBRA modeling

Modified yeast 7.0 modeling with updated biomass function

Under the consideration that a modeling tool with a more comprehensive biomass function could provide more information, we repeated modeling of the 5 relevant enzymes using a biomass function with added requirements for ferroheme b and charged tRNA(Tyr) in the mitochondrion. This model was able to implicate three additional proteins with function in essential reactions: delta-aminolevulinic acid dehydratase (ALADH), P05373; mitochondrial ferrochelatase, P16622; and mitochondrial tyrosine-tRNA ligase, P48527. This additional examination, in conjunction with sequence and pathway conservation, supports our argument for inclusion of the pertinent pathways in yeast cell models. A summary of the findings can be seen in Table 5. (see Additional file 1: Table A3 for full output).

The above Yeast 7.0 results, combined with literature support, sequence and pathway annotations, maintain our hypothesis that we should be able to use yeast orthologs with active site conservation to model the effect of a variation at the active site of its human ortholog in order to obtain a predictive overview of altered enzyme activity resulting from the variation. Furthermore, we see the potential for the utility of this method to improve as advanced modeling efforts continue to progress.

Biomedical relevance

We hypothesize that nonsynonymous variation in an enzyme’s active site should have deleterious effects on the enzyme’s activity and function. Although only one of the six human enzyme SNPs has a direct link to the Online Mendelian Inheritance in Man® database (OMIM®) [41] in dbSNP (rs28936396 in P48637, human glutathione synthetase linked to OMIM 266130 – Glutathione synthetase deficiency), all six have literature related to dysfunction and disease caused by variants found at other positions within the protein sequence. If our assumptions regarding active site variation are correct, then nsSNVs at the active site of enzymes involved in metabolic syndromes will likely alter protein function and contribute to similar syndromes. Furthermore, yeast orthologs for these enzymes with strong pathway and active site conservation may be used to verify the metabolic effects of such variation with respect to metabolite concentrations. Although the current Yeast 7.0 model supports this hypothesis for RPIA alone, we argue for inclusion of the pathways and metabolites associated with these proteins in future yeast models based on consideration of the augmented model findings, literature and conservation among orthologs. Thus, we present five potential disease-related proteins with active site nsSNVs (along with yeast orthologs) in need of laboratory experimentation to validate the veracity of disease related to the SNV.

Human P49247/Yeast Q12189

Ribose-5 phosphate isomerase (RPI) participates in the pentose phosphate pathway (PPP), an alternative pathway for glucose oxidation accounting for up to 20% of glucose oxidation in normal tissue. [42] Synthesis of ribose 5-phosphate from ribulose 5-phosphate via RPI is required for nucleotide and nucleic acid synthesis, and participates in the downstream production of glycolytic intermediates [43]. With reduced RPI activity, ribose 5-phosphate becomes unavailable and the necessary glycolytic metabolites are not produced. Enzyme deficiency results in clinical symptoms of leukoencephalopathy and peripheral polyneuropathy. The RPIA gene has also been identified to be hypermethylated at CpG sites in breast cancer [44].

As of 2004, RPI was only the second known inborn error in the reversible portion of the PPP [45]; the PPP connects pentose phosphates to glycolytic intermediates. Furthermore, there are currently only 2 mutations listed in the Human Gene Mutation Database [46] (accessed 8/22/13) for the RPIA gene, neither of which is the aspartic acid to tyrosine variation at active site position 160 reported here (rs11549730). RPI has been observed to be variably conserved among species ranging from 52% similarity between human and E. coli to strong conservation noted between human and other mammals [47] and RKI1 disruption has been shown to be deleterious to yeast growth [48].

Human P13716/Yeast P05373

Delta-aminolevulinic acid dehydratase (ALADH), encoded by the ALAD gene in humans, catalyzes the first common step in heme and other tetrapyrrole biosynthesis [49] and has also been observed to interact with the proteasome [50]. The wild-type enzyme exists as a high-activity homooctamer, but alternative lower-activity hexamer assemblies have been found [51]. In fact, naturally occurring ALAD porphyria-associated human variants have been shown to have an increased susceptibility to hexamer-stabilizing inhibitors. Lead-poisoning is also hallmarked by reduced ALAD activity. 5-aminolevulinate accumulates as a result of the lowered activity and is responsible for the toxic effects observed in both disease states. It has been suggested that small molecule stabilization of alternate oligomers (morpheeins) can be used as a therapy to target delta-aminolevulinic acid dehydratase [52], in addition to human immunodeficiency virus integrase, tumor necrosis factor α, mammalian ribonucleotide reductase and other disease-related proteins for which there is evidence of alternate, functionally differential oligomers [53].

Mutants in the corresponding yeast HEM2 gene have been found to disrupt heme biosynthesis, requiring provision of heme for growth in culture [54]. Yeast and human ALADH differ with respect to Zn cofactor binding, thus care is required to ensure homology models account for this difference in stoichiometry [55]. However, comparison of primary sequences across a diverse array of species shows a strong degree of similarity in the region of the active site and a zinc-binding motif [56].

There are many studies discussing ALAD-porphyria associated mutations [57, 58] but responsibility for the disease state has not yet been assigned to the deletion at nucleotide position 843. Interestingly, this deletion in the active site codon also results in a frameshift, further complicating comprehension of potential disease mechanism(s).

Human P22830/Yeast P16622

Encoded by the FECH gene, ferrochelatase is the final enzyme in heme biosynthesis, catalyzing the formation of protoheme IX from protoporphyrin IX and ferrous iron [59]. Ferrocheletase is important to human health due to the critical and diverse functions of its product cofactor, heme. Defects in this enzyme have been related to erythropoietic protoporphyria (EPP). EPP is a relatively benign disease with predominant manifestation of photosensitivity, although it has been reported with liver complications in 2-5% of patients [60]. Biological findings include elevated concentrations of protoporphyrin in erythrocytes, plasma and feces, and a large amount of protoporphyrin in the skin [61].

Comparisons between ferrochelatases of multiple species reveal the proteins belonging to human, mouse, chicken, frog and Drosophila melanogaster to be metalloenzymes with a [2Fe-2S] cofactor, thought to play a structural stabilizing role in humans [62]. The corresponding bacteria, yeast and plant proteins contain no iron-sulfur center, but the ferrochelatase from S. cerevisiae is similar in other respects to eukaryotic ferrochelatases. Yeast mutants with total or partial FECH activity deficiency and resultant protoporphyrin accumulation have been isolated [63].

Human Q9Y2Z4/Yeast P48527

Aminoacyl-tRNA synthetases (ARSs) catalyze the linkage of specific amino acids to their cognate tRNAs [64]. Mitochondrial tyrosyl-tRNA synthetase (TyrRS) is responsible for charging tyrosine to its tRNA. Although belonging to class I synthetases, tRNA recognition occurs in a way more common to class II ARSs [65]. Defects in TyrRS are associated with myopathy with lactic acidosis and sideroblastic anemia 2 (MLASA2). It has been proposed that diminished aminoacylation activity of the enzyme leads to decreased mitochondrial protein synthesis and mitochondrial respiratory chain dysfunction [66]. TyrRS also has known involvement with cell-signaling [67] and angiogenesis [68]. Although there are multiple known disease-related TyrRS mutants, it is unclear for some mutants whether pathogenicity results from dysfunction in aminoacylation or cytokine activation.

Alignments with TyrRSs from several species reveals 34% identity with yeast and similarly low identity across the board, with greatest conservation confined to the N-terminus [69]. Gly-244 and Asp-246 are conserved among class I synthetase catalytically important residues, but a total of 8 residues are known to participate in hydrogen bonding [70].

Human Q96GX9/Yeast P47095

Methylthioribulose-1-phosphate dehydratase, encoded by the APIP (APAF1-interacting protein) gene, catalyzes the dehydratase step of the methionine salvage pathway, converting methylthioribulose-1-phosphate into 2,3 dioxomethylthiopentane-1-phosphate [71]. The methionine salvage pathway recycles sulfur metabolites by recovery of methionine from methylthioadenosine (MTA). This pathway plays a critical role in many biological functions including apoptosis, inflammation and cancer [72].

APIP is highly expressed in skeletal muscle and known for an anti-apoptotic role [73]. One common APIP SNP was observed to reduce enzyme activity, resulting in an increase in cell death in response to Salmonella[74].

Currently excluded from the Yeast 7.0 model, the MDH1 enzyme is of critical importance to the pathway. Loss of function leaves the cell with no isozymes for methionine salvage, as demonstrated in the YeastCyc pathway/genome database (PGDB) within BioCyc [75]. The enzyme and the entire methionine salvage pathway seem to be functionally conserved in yeast [76] and should therefore serve as an adequate model in yeast experimentation to validate predicted loss of function due to the contained nsSNV.

Datasets and annotated information

Proteomes and all protein subsets were retrieved from UniProtKB/Swiss-Prot (UniProt release 2012_06 [28]. nsSNV data were extracted from UniProtKB/Swiss-Prot, dbSNP (build 135) [29], COSMIC (version 59) [77] and NCI 60 [78] (accessed June 28, 2012). The SNVDis tool (https://hive.biochemistry.gwu.edu/snpdis.cgi?cmd=dmSnpdis) [12] was used to retrieve the set of all proteins with nsSNVs at the active site (as annotated by UniProtKB/Swiss-Prot or CDD). Pathway information was retrieved from Kyoto Encyclopedia of Genes and Genomes (KEGG) [79] using the Protein Analysis Through Evolutionary Relationships (PANTHER) Classification System [80]. Substrate, product and disease information was gathered manually from KEGG and UniProtKB/Swiss-Prot. Official names and PubChem [81] Compound IDs (CIDs) for compounds with specific chemical structures, or Substance IDs (SIDs) for those without, were recorded for all substrates and products.

Analysis

Reviewer 1: Dr. Daniel Haft

Previous work by this group identified 559 human enzymes in which missense mutations replace an enzyme’s active site residue in rare alleles. Of these, 34 are assigned a single enzymatic activity rather than a panel of alternate substrates, with the reaction belonging to a single pathway in human as documented by KEGG. For this paper, the collection is winnowed further, to count only those for which bi-directional best BLAST matches between human and yeast identify putative ortholog pairs, uncomplicated by paralogs that could overlap in function. This stringent filter gets the list down to six enzymes, the topic of this paper. A tool for flux balance analysis (FBA) is available for yeast, but not for human. Thus, for these six enzymes, computation based on a yeast metabolic model leads to predictions that loss of the same function in human causes disease.

In one sense, there is no surprise here. Mutations away from the active site already were known to cause human disease for five of the six enzymes. Loss of the active site residue should be no less catastrophic to function than truncation or frameshift errors. It is possible that a heterozygous active site mutant could produce a currently unrecognized form of a disease now known only for homozygotes of an incompletely inactivated enzyme, but that is not the point of the paper. The yeast model is haploid. Presumably, the implications are for homozygotes of the mutant form in human.

Authors’ response: Perhaps it wasn’t clear from our original discussion, but we agree that the haploid yeast model is only suited to truly “model” homozygous cases. However, we can envision potential scenarios in which a haploid model could still provide predictive quantitative value in assays in highly conserved pathways for loss of function due to heterozygous mutation at active sites in humans. We did revisit the text relevant to this issue and modified it to improve clarity.

Furthermore, we are not claiming that homozygous alleles are altogether too rare, in fact the journal article we cite here (Ng PC1, Levy S, Huang J, Stockwell TB, Walenz BP, Li K, Axelrod N, Busam DA, Strausberg RL, Venter JC. Genetic variation in an individual human exome. PLoS Genet. 2008 Aug 15;4(8):e1000160) reports that just under 50% of all nonsynonymous SNPs in HuRef are homozygous. However, they go on report that only about 15% are likely to affect protein function in a deleterious way and in this case, heterozygous nsSNPS are two times more likely than homozygotes to be predicted to affect protein activity. Similarly, rare nsSNPS are twice as likely as common nsSNPs to affect protein function. Thus, at the active site and other functionally important residues, protein-affecting SNPs should be selected against and should be less likely to be observed in homozygote form.

At your suggestion, we conducted a new PubMed search and were able to find an article identifying a homozygous mutation in the catalytic domain of the MMP2 gene in two sisters with Winchester syndrome (Rouzier C1, Vanatka R, Bannwarth S, Philip N, Coussement A, Paquis-Flucklinger V, Lambert JC. A novel homozygous MMP2 mutation in a family with Winchester syndrome. Clin Genet. 2006 Mar;69(3):271–6.) The long-term combination of (1) improved classification and annotation of variation data and (2) simple increase in available studies should allow us to find these cases and truly study the quantitative differences of heterozygous and homozygous states, ultimately helping to develop and refine possible detection assays.

The surprising finding, actually, seems to be that even under the stringent filtering used – conserved orthologs from yeast to human of enzymes in a conserved pathway – five of the six genes nominated for study are associated with described human genetic diseases, not simply non-viability. The alleles responsible for the genetic disorders seen in patients may encode enzymes that retain some residual function, probably more than an active site mutant retains, allowing for viability.

Rather than presenting a new method to find previously unrecognized human genetic disease, this paper forms a field test of the bioinformatics infrastructure one would need for such a pipeline in the future. If FBA is not yet available for human cells, the yeast model can stand in for a few validated human/yeast orthologs. When it turns out that zeroing activity of several genes causes no change, it suggests the FBA model needs repair. The authors describe making such repairs, on their way to completing their proof-of-principle for notion that a yeast metabolic model, and an infrastructure for data integration, could point the way to the prediction of diseases that could be caused by certain missense mutations.

Authors’ response: Thank you very much for your detailed comments and review which helped us to focus our editing efforts and improve our arguments.

Reviewer 2: Dr. Igor B. Rogozin

The authors analyzed nonsynonymous single nucleotide variations (nsSNVs) in active site residues in human-yeast protein alignments. This is a promising approach. I do not see major methodological problems. The only problem that I see is a discussion about "loss of function via haploinsufficiency or dominant negative phenotypes". I understand that this is a complex problem. For example, the vast majority of disease-associated mutations in the RPE65 gene (Poliakov E, Gubin AN, Stearn O, Li Y, Campos MM, Gentleman S, Rogozin IB, Redmond TM. Origin and evolution of retinoid isomerization machinery in vertebrate visual cycle: hint from jawless vertebrates. PLoS One. 2012;7: e49975) are recessive mutations. However a dominant mutation was recently found in this gene. This and many other examples (including those described in the paper) suggest that the problem of dominant mutations is extremely complex. I am not sure that the fraction of such mutations is high although I am not aware of such estimate in active site residues.

Authors’ response: Thank you for your positive comments regarding our approach and your general suggestions as well. We do realize the complexity of this issue and the necessity of exercising caution when applying models across organisms without knowing the full characterization of the deleterious outcome. As I mentioned in my response to Dr. Haft above, we have discussed the issue of zygosity, phenotype and loss-of-function both prior to the initial submission and since receipt of your responses. We do not want to ignore the modeling ramifications of differential ploidy, but we also don’t want to exclude candidate proteins from downstream quantification steps as it may be possible that future in vivo assays (developed as a product of the broader study) can still have some impact on heterozygous disease states even if in silico models cannot properly simulate them. We have carefully reread and edited these sections in response to your comments and we feel that we have clarified our reasons for our treatments of these issues.

Reviewer 3: Dr. I. King Jordan

This reviewer provided no comments for publication.