Substrate Channeling via a Transient Protein-Protein Complex: The case of D-Glyceraldehyde-3-Phosphate Dehydrogenase and L-Lactate Dehydrogenase

Substrate channeling studies have frequently failed to provide conclusive results due to poor understanding of this subtle phenomenon. We analyzed the mechanism of NADH-channeling from D-glyceraldehyde-3-phosphate dehydrogenase (GAPDH) to L-lactate Dehydrogenase (LDH) using enzymes from different cells. Enzyme kinetics studies showed that LDH activity with free NADH and GAPDH-NADH complex always take place in parallel. The channeling is observed only in assays that mimic cytosolic conditions where free NADH concentration is negligible and the GAPDH-NADH complex is dominant. Molecular dynamics and protein-protein interaction studies showed that LDH and GAPDH can form a leaky channeling complex only at the limiting NADH concentrations. Surface calculations showed that positive electric field between the NAD(H) binding sites on LDH and GAPDH tetramers can merge in the LDH-GAPDH complex. NAD(H)-channeling within the LDH-GAPDH complex can be an extension of NAD(H)-channeling within each tetramer. In the case of a transient LDH-(GAPDH-NADH) complex, the relative contribution from the channeled and the diffusive paths depends on the overlap between the off-rates for the LDH-(GAPDH-NADH) complex and the GAPDH-NADH complex. Molecular evolution or metabolic engineering protocols can exploit substrate channeling for metabolic flux control by fine-tuning substrate-binding affinity for the key enzymes in the competing reaction paths.

www.nature.com/scientificreports www.nature.com/scientificreports/ NADH channeling can be one of the key functions in the supramolecular organization between glycolytic enzymes. NADH channeling has been explored with limited success in the last 40 years 7,15,[18][19][20][21][22] . A network of channeling reactions between different NAD(H) dehydrogenases can maintain the separation between different NAD+/NADH pools in cells 1,8,13,23,24 . Such separations can be the key mechanism in the regulation of cell physiology and tumor development 1,8,13,23,24 . Cytosolic and mitochondrial NAD+/NADH pools are highly un-equilibrated 25 , and still closely integrated in control of cellular energy production 8,26 . Cytosolic NAD+/NADH pool can also affect the cytosolic NADP+/NADPH pool through dual-specificity malate dehydrogenase 27 . The separation between NADH and NADPH pools can regulate the balance between catabolic and anabolic processes in cells 8,26 .
GAPDH is the most abundant cytosolic dehydrogenase, that binds the majority of cytosolic NAD(H) molecules 1,2,13,27,28 . NADH channeling from GAPDH to different NADH-dehydrogenases could regulate the separation between different NAD(H) pools in cells 8,26 . The changes in NAD+/NADH concentrations have a peculiarly strong effect on GAPDH structure 1,2,13,[27][28][29][30] . A decrease in NADH concentration leads to partial dissociation of the GAPDH tetramers 2,13,31 . NAD(H) binding to different subunits in GAPDH tetramer results in allosteric regulation of GAPDH activity and strong negative cooperativity that can change the NAD(H) binding affinity by thousand-fold [28][29][30] . These effects are species-specific and known for decades, but their physiological significance is still unknown 2,13,31 .
LDH-GAPDH complex was observed in cell extracts 20 , and with purified proteins in conditions that mimic high protein concentrations in cytosol 21 . Such LDH-GAPDH complex can regulate the major metabolic pathways, however, functional consequences of LDH-GAPDH interaction have not been investigated. Here we have provided a new set of evidence that shows that LDH and GAPDH tetramers can form a transient supramolecular complex that can simultaneously support channeled and diffusive reactions. Such competition between channeled and diffusive paths can fundamentally change our understanding of metabolic regulation and metabolic flux control within a glycosome. Most notably, changes in the substrate-binding affinity can be a molecular mechanism for fine-tuning of metabolic pathways in enzyme evolution and metabolic engineering protocols. Our results also show that it is impossible to design experiments that can conclusively analyze substrate channeling in cells if we do not understand the underlying molecular principles and the properties of the related enzymes.

Results
Molecular dynamics simulations of the interaction between rabbit muscle GAPDH and rabbit muscle LDH. LDH-GAPDH complex can be observed in cell extracts 20 , and with purified proteins in conditions that mimic high protein concentrations in cytosol 21 . We used different structure analysis techniques to calculate the structure of rmLDH-rmGAPDH complex. Distribution of electric potentials on the protein surface can give initial insights into the structural elements that can support substrate channeling [32][33][34] . We find that NADH binding sites on two adjacent monomers in tetramers of rmLDH 35 are connected with a 19.2 Å long groove that is dominated by the positive potential (Fig. 1A,B). The groove is observed only when the active site loops on rmLDH are open (Fig. 1A, residues 98-105 35 ). In rmGAPDH tetramers 36 , the NADH binding sites on two adjacent monomers are enclosed within 18.4 Å nm wide gulf between the monomers (Fig. 1C). The gulf is entirely www.nature.com/scientificreports www.nature.com/scientificreports/ dominated by positive electric fields (Fig. 1C,D). The electric fields in NAD(H) binding sites are projecting dominant positive electric fields in the space around each protein (Fig. 1B,D). Such positive fields can limit diffusion of negatively charged NAD(H) molecules between the adjacent subunits and facilitate substrate channeling within rmLDH or rmGAPDH tetramers 32,33 . The positive fields around each protein can visibly decrease when the proteins bind negatively charged NAD(H) molecules ( Fig. 1B-D).
Based on the calculated surface potentials, we have presented a possible binding orientation for the rmLDH-rmGAPDH complex in which four NAD(H) binding sites are enclosed in one positive field ( Fig. 2A). The number of possible docking orientations is limited by the two main requirements 32,33 . The distance between the channeling sites should be between 2 to 5 nanometers, and the channeling path should be continuous and at least partially secluded from the surrounding solution.
Multiscale MD calculations and molecular docking studies showed that rmLDH and rmGAPDH can form a dynamic complex facing each other with their NAD(H) binding sites (Figs. 2 and 3, Supp. Figures 1-8 and 10, and Supp. Video 1). The complex breaks apart when the two enzymes are saturated with NAD(H) molecules (Supp. Video 1). This is consistent with the earlier experimental studies of the interaction between LDH and GAPDH molecules which showed that saturation with NADH leads to complex breakdown 21 . When the enzymes are saturated with NAD(H) they form some random contacts but ultimately slip apart (Supp. Video 1, and Supp. Fig. 8G-I). Superimposed structures of NADH-free and NADH-bound forms of each enzyme showed that NAD(H)-bound structures cannot support the formation of the complex due to the repositioning of the interacting amino acids (Supp. Video 2). NADH binding to rmLDH will cause the closing of the active site cleft 35,37 . NAD(H) binding to GAPDH will cause 7.4 Å contraction between the two subunits that interact with LDH 36 .
Repeated MD simulations showed that the conformational changes around NAD(H) binding sites lead to stochastic differences in the rate of interaction build-up and in variability in the number of binding interactions (Supp. Figs. 7, 8 and 10). The structures from the MD frames that show differences in the number of binding interactions have been superimposed to analyze the interaction mechanism. The highest number of binding interactions can be observed when both rmLDH and rmGAPDH tetramers are present in their open conformations that are dominant in NADH-free structures ( Fig. 3 and Supp. Videos 1 and 2). Most notably stochastic differences in the repeated simulations can be attributed to the random wobbling in the position of the active site loop on rmLDH (a.a. 98-104). When the active site loop is closed the interacting proteins start wobbling around the axis that is perpendicular to the plane of interaction (Supp. Videos 1 and 2). The delicate nature of the presented LDH-GAPDH interactions can explain why it was so difficult to measure that complex 20,21 . The LDH-GAPDH tetramers form a symmetric complex that can support the association of multiple LDH and GAPDH molecules in a polymer, which can explain the poor solubility of the complex 20,21 .
When rmLDH and rmGAPDH form a complex, the positive cavities on the surface of each enzyme merge to form a central positive cavity under the protein surface ( Fig. 2C-E). The cavity connects four NAD(H) binding sites with an average separation of 2.9 nm between the adjacent sites ( Fig. 2C-E). Thus, NADH channeling within the rmLDH-rmGAPDH complex can be an extension of NADH channeling between the two adjacent monomers in rmLDH and rmGAPDH tetramers ( Figs. 1 and 2). www.nature.com/scientificreports www.nature.com/scientificreports/

Interaction between LDH and GAPDH molecules from different cells. LDH and GAPDH
molecules from different cells have highly conserved structures but differ in the net charge and the NAD(H) binding mechanism 13,21,28,38 . Both rabbit muscle and porcine heart LDH tetramers can form complex with rmGAPDH 21 . phLDH has 75% sequence identity and 93.1% sequence similarity with rmLDH (5LDH 39 and 3H3F 35 ). The two molecules have opposite net charge (pI(rmLDH) = 8.1, pI(phLDH) = 4.6) 21 ). The corresponding structures can overlap with the RMSD value of 1.73 Å. The amino acids that form binding interactions in the rmLDH-rmGAPDH complex appear to be conserved between rabbit muscle LDH and porcine heart LDH (Fig. 3). Repeated MD simulations showed that phLDH-rmGAPDH complex can form in average 2 more binding interactions than the rmLDH-rmGAPDH complex (Supp. Fig. 9A-C). The most notable difference is the replacement of Gln 228 in rmLDH with Glu 227 in phLDH (Fig. 3). The substitution is on a flexible helix that can adapt to different orientations between the interacting proteins ( Fig. 3). Porcine and rabbit heart LDH molecules share 99.1% sequence similarity, 95.8% sequence identity, and 100% identity at the interaction sites (Fig. 3).
Rabbit muscle GAPDH and baker's yeast GAPDH have 65.5% sequence identity and 84.8% similarity (PDB: 3PYM, isozyme 2 36,40 ). The corresponding structures can overlap with the RMSD value of 0.556 Å. byGAPDH and rmGAPDH have very different NADH binding affinities and opposite net charge (pI(byGAPDH) = 6.5, pI(rmGAPDH) = 8.2) 36,40 . Repeated MD calculations showed that byGAPDH-rmLDH complex can form in average 8 ± 2 binding interactions, i.e., two interactions less than rmGAPDH:rmLDH complex (Fig. 3, Supp. Enzyme buffering tests with different LDH molecules as enzyme acceptor and GAPDH molecules as NADH donor enzyme. We have designed enzyme assays that can measure NADH channeling from GAPDH-NADH complex to different LDH isozymes (Supp. Figs. [11][12][13][14]. The measurements were designed to mimic LDH activity in the cytosol, where the majority of NAD(H) molecules are bound to GAPDH, the most abundant dehydrogenase in cells with the highest NADH binding affinity 13,27 . The enzyme assays with two enzymes that share a common substrate can be challenging to design and interpret 18,19 , especially in the case of transient LDH-GAPDH interaction (Supp. Figs. 11 20,21 ) Thus, we used extensive numerical simulations 41 in preparation for optimal assay design and data interpretation (Supp. Figs. 11 to 14). The simulations showed that the free-diffusion and the channeling paths always take place in parallel. Thus, we measured LDH activity using two approaches that can modulate the relative contributions from the two paths: increasing GAPDH concentration with fixed NADH concentration and decreasing NADH concentration with fixed GAPDH concentration. In the case of no channeling the measured activity will be equal to the calculated LDH activity on free NADH substrate (Eqs. 1-2, Supp. Fig. 11). The expected free-diffusion activity can be calculated using Michaelis-Menten constants for LDH activity with its NADH substrate (Table 1) and dissociation constant Kd for GAPDH-NADH complex ( Table 2). In the case of channeling, the measured LDH activity will be higher than the calculated activity on the free NADH substrate (Supp. Fig. 11). We found with four different enzyme pairs that LDH activity in the presence of a large excess of GAPDH is significantly higher than the expected activity on free NADH substrate (Figs. 4 and 5).
In the first set of experiments, we compared rmLDH and phLDH as enzyme acceptors with rmGAPDH as a donor enzyme (Fig. 4A,B). LDH activities were measured with fixed NADH concentration (40 μM) and variable rmGAPDH concentrations (100 to 240 μM of NADH binding sites or 3.5-8.5 mg/ml protein) (Fig. 4A). In these conditions, more than 99% of all NADH molecules are bound to rmGAPDH (Supplement Table 2). The free www.nature.com/scientificreports www.nature.com/scientificreports/ NADH concentration is at least 11-fold lower than the related K M constants for each LDH molecule (Supplement Table 2). Increase in rmGAPDH concentration leads to an increase in the ratio between measured and calculated activity from 1.2 to 1.85 fold (eqn. 3 in methods). In the same experiments with phLDH, increase in rmGAPDH leads to an increase in the ratio between the measured and calculated activity from 1.3 to 2.0 fold (eqn. 3 in methods). The experiments with fixed rmGAPDH concentration showed that NADH decrease from 40 to 10 μM leads to increase in the ratio between measured and calculated activity from 1.65 to 2.25 fold for rmLDH and 1.7 to 2.28 fold for phLDH ( Fig. 4B and Supplement Table 2).
In the second set of experiments, we have replaced rmGAPDH as enzyme donor with byGAPDH ( Fig. 5A,B). The activities or rmLDH and phLDH were measured at fixed NADH concentration (40 μM) and variable concentrations of byGAPDH (240 to 480 μM of NADH binding sites or 8.5-17.6 mg/ml protein). In those conditions, more than 96% of NADH molecules are bound to byGAPDH (Supplement Table 2B). The concentration of free NADH is at least five-fold lower than the respective Km constants for rmLDH or phLDH (supplement Table 2B). In experiments with rmLDH, increase in byGAPDH concentration leads to an increase in the ratio between measured and calculated activity from 3.5 to 5.8 fold (Fig. 5A and Supplement Table 2B). In the same experiments with phLDH, increase in byGAPDH concentration leads to an increase in the ratio between measured and calculated activity 3.5 to 6.1 fold ( Fig. 5A and Supp. Table 2B). The experiments with fixed byGAPDH concentration showed that NADH decrease from 40 to 10 μM leads to increase in the ratio between measured and calculated activity from 4.85 to 6.2 fold for rmLDH and 5.2 to 6.6 fold for phLDH (Fig. 5B Supplement Table 2).
In conclusion, consistent with the molecular dynamic studies, all four experiments showed that the channeled and the diffusive reactions always take place in parallel. The channeling can be measured only at extremely high protein concentrations which mimic cytosolic conditions that favor the formation of GAPDH-NADH complex at the expense of free NADH (Figs. 4 and 5) 1,2,13,27,28 . These results are consistent with the earlier experiments 20,21 and with the numerical simulations (Supp. Figs. [11][12][13][14]. We find that the highest GAPDH concentrations are limited by the enzyme solubility, while the lowest NADH concentrations were limited by assay sensitivity (Supplement Table 2). In all four experiments the conditions that favor the formation of GAPDH-NADH complex lead to a decrease in measured LDH activity, about 38% decrease with rmGAPDH and about 10% decrease with byGAPDH (Figs. 4 and 5). The decrease can be attributed to the shift from the faster turnover in free-diffusion reaction to the slower turnover in channeling reaction (Supp. Figs. 11 to 14). The turnover rates in channeling reactions are always slower than the turnover rates for the diffusive reaction due to one extra step: NADH dissociation from GAPDH within (GAPDH-NADH)-LDH complex (Supp. Figs. 12 and 13). In all four measurements, the apparent K M values for LDH molecules with GAPDH-NADH substrate are higher than the K M values for free NADH (Tables 1,3, and Figs. 4B and 5B). Thus, in the case of channeling the extent of saturation of LDH activity with the NADH substrate is a result of competition between the channeled and diffusive reaction paths (Supp. Fig. 14).
Control enzyme buffering experiments: competition assays with anti-phLDH antibodies. In a negative control experiment, the enzyme buffering experiments were repeated in the presence of an excess of polyclonal anti-LDH antibodies (Fig. 6). The large anti-LDH IgG molecules have little effect on the binding of small NADH molecules to phLDH. On the other hand, approximately 45% inhibition was observed in the presence of rmGAPDH (240 μM NADH binding sites or 8.5 mg/ml protein). Similarly, about 80% inhibition was observed in the presence of byGAPDH (450 μM NADH binding sites or 15.9 mg/ml). The anti-LDH antibodies have a much higher binding affinity for LDH molecules than the GAPDH molecules, and thus the anti-LDH antibodies can readily interfere with the formation of different LDH-GAPDH complexes (Fig. 6). www.nature.com/scientificreports www.nature.com/scientificreports/ Figure 4. (A,B) The activity of rabbit muscle LDH and porcine heart LDH was measured in the presence of a large excess of rabbit muscle GAPDH. (A) steady-state activities of rmLDH (10 nM) or phLDH (17 nM) were measured at fixed NADH concentration (40 µM) in the presence of increasing concentration of rmGAPDH (100 to 240 μM of NADH binding sites). Increase in rmGAPDH concentration leads to the disproportional decrease in the measured and the calculated free-diffusion activities (lower panels) which results in the increase in the ratio between the two activities (upper panels, and Eq. 2−3). These results indicate that LDH molecules can use GAPDH-NADH complex as a substrate in addition to free NADH, i.e. NADH channeling (Supp. Fig. 11). (B) steady-state activities of rmLDH (10 nM) or phLDH (17 nM) were measured in the presence of decreasing NADH concentration with rmGAPDH fixed at 200 μM (NADH binding sites). The decrease in NADH concentration leads to the disproportional decrease in the measured and the calculated free-diffusion activities (lower panels), which results in increase in the ratio between the two activities (upper panels, and Eq. 2-3). The red curve represents the Michaelis-Menten profile for LDH activity with the rmGAPDH-NADH complex as the substrate, which was calculated by subtracting the calculated free-diffusion profile from the measured profiles (Supp. Fig. 14). The calculated apparent K M constant for rmLDH is 35 ± 5.5 μM and 59 ± 6 μM for phLDH (Table 3). Thus, the observed K M constants are a result of competition between the channeled and free-diffusion paths (Supp. Fig. 11). . Increase in byGAPDH concentration leads to the disproportional decrease in measured and calculated free-diffusion activities (lower panels) which results in increase in the ratio between the two activities (upper panels, and Eq. 2-3). These results indicate NADH channeling from byGAPDH-NADH complex to rmLDH or phLDH (Supp. Fig. 11). (B) steady-state activities of rmLDH (10 nM) or phLDH (17 nM) were measured in the presence of decreasing NADH concentration with byGAPDH fixed at 480 μM (NADH binding sites). The decrease in NADH concentration leads to the disproportional decrease in the measured and the calculated free-diffusion activities (lower panels), what results in the increase in the ratio between the two activities (upper panels, and Eq. 2-3). The red curve represents the Michaelis-Menten profile for LDH activity with the byGAPDH-NADH complex as the substrate, which was calculated by subtracting the calculated free-diffusion profile from the measured profiles (Supp. Fig. 14). The calculated apparent K M constant for rmLDH is 78 ± 3 μM and 116 ± 10 μM for phLDH (Table 3). Thus, the observed K M constants are a result of competition between the channeling and free-diffusion paths.
www.nature.com/scientificreports www.nature.com/scientificreports/ At the highest antibody concentration, the measured phLDH activity in the presence of GAPDH molecules is close to calculated activity for the reaction that depends on free diffusion (methods Eqs. 1-2). A good agreement between the calculated and the measured activity indicates that the K D constant for GAPDH-NADH complex and Michaelis-Menten parameters for LDH activity were measured with high accuracy (Tables 1 and 2).
These results support our proposals that the channeled and the diffusive reactions always take place in parallel.

Analysis of interaction between phLDH and byGAPDH by analytical ultracentrifugation.
Our earlier studies of the interaction between phLDH and rmGAPDH showed that AUC studies can be a very sensitive method for the detection of interaction between LDH and GAPDH molecules 21 . The same approach is now used for the analysis of the interaction between phLDH and byGAPDH, the two enzymes that show channeling (Figs. 4 and 5). The sedimentation constant for phLDH is s 0 20,w = 7.5 ± 0.2 S, while apparent molecular mass is Mr = 143,6 ± 0.8 kDa. The calculated Mr is within 3% of the calculated Mr based on amino acid sequence (146,5 kDa). For byGAPADH the sedimentation constant s 0 20,w = 7.65 ± 0.2 S, while apparent molecular mass is 142.1 ± 0.7 kDa what is to within 1.4% of the Mr value that can be calculated based on amino acid sequence (143 kDa).
The sedimentation profiles for phLDH and rmGAPDH have such high similarity that that two proteins mixture can give a good fit to one component model, with the best fit values s 0 20,w = 7.6 ± 0.4 S, Mr = 143,6 ± 0.9 kDa (Fig. 7A). The good fit to the one-component model shows that there is no detectable complex between phLDH-byGAPDH in the two-protein mixture. A good fit to one component model can be also used to estimate the sensitivity of our AUC experiments to detect traces of association. The calculated s and Mr values for a single component model can be used in Claverie simulations in SEDFIT program to simulate sedimentation profile for the self-dimerization model assuming 10% association 42 . In our experiment, each protein in the mix is present in the total concentration of 6 μM, thus if there is 10% association the complex concentration will be 0.6 μM, the concentration of free protein is 5.4 μM and the corresponding dissociation constant K D is 48.6 μM (i.e. K D = (5.4)*(5.4)/0.6). The measured and the simulated profiles show detectable differences (Fig. 7). We conclude that there is no interaction between phLDH and byGAPDH with dissociation constant lower than 48.6 μM.   www.nature.com/scientificreports www.nature.com/scientificreports/ Presented results (Figs. 4-7) and the earlier LDH-GAPDH interaction studies 20,21 indicate that NADH channeling does not require a high affinity complex between GAPDH and LDH. These results indicate that the channeling is taking place within a transient LDH-(GAPDH-NADH) complex (Supp. Fig. 11 and 43 ).

Discussion and conclusions
The molecular mechanism and the regulation of NADH channeling from GAPDH-NADH complex. Mechanisms that regulate the integration of different enzymes in cellular metabolic pathways are among some of the most fundamental unanswered problems in biochemistry. We have observed differences between rmGAPDH and byGAPDH that can give some general insights into the channeling mechanism and the enzyme buffering experiments (Figs. 4-5, Supp. Tables 2 and 3). rmGAPDH and byGAPDH have very similar structures but differ in their NAD(H) binding mechanism 13,29,40,44 . Briefly, rmGAPDH (PDB: 1J0X) and byGAPDH (PDB: 3PYM, isozyme 1) have 65.5% sequence identity and 84.8% similarity. The corresponding structures can overlap with RMSD value of 0.556 Å. The two enzymes have such high structural similarities that they can even form heterotetramers 13,29,40,44 . The heterotetramers readily fall apart after saturation with NAD(H) due to different conformational changes caused by the NADH binding to yeast and rabbit subunits 13,29,40,44 . rmGAPDH shows strong negative cooperativity in NADH binding to different subunits in the tetramer 13,44 . byGAPDH does not show negative cooperativity in NADH binding 40 and its binding affinity is 10 times weaker than the NADH binding affinity for rmGAPDH (Table 1).
Numerical analysis showed that in the case of transient LDH-(NADH-GAPDH) interaction weaker NADH binding to GAPDH can favor channeling (Supp. Figs. 11-13). The calculations with increasing off-rates showed that in the case of transient protein-protein interactions the channeling depends on an overlap in timing of two events: dissociation rate for LDH-(NADH-GAPDH) complex and NADH dissociation rate for GAPDH-NADH complex (Supp. Figs. 11-13). The higher K D constant for the byGAPDH-NADH complex is in large part a result of higher off-rates for GAPDH-NADH complex 29,[44][45][46][47][48] . Numerical analysis showed that the increase in off-rates for GAPDH-NADH complex can qualitatively reproduce experimentally observed differences between rmGAPDH and byGAPDH molecules (Supp. Figs. 12 and 13).
The correlation between NADH channeling and dissociation constant for GAPDH-NADH complex can indicate that the negative cooperativity in NADH binding to different subunits in GAPDH tetramer could regulate the NADH channeling 1,2,44 . Allosteric regulation of substrate channeling has been reported in some of the earlier studies 17,33,49 . The negative cooperativity in NADH binding to different subunits in GAPDH molecules is, in essence, a type of allosteric regulation 29,46,50 . The negative cooperativity can affect NADH binding affinity to different subunits in GAPDH tetramer by about 1000-fold 47,48 , from 10 nM to 30 μM. The 1000-fold difference in the binding affinity can be in large part a result of differences in the off-rates for the GAPDH-NADH complex 1,44,45 . In other words, we hypothesize that we can measure no-channeling, low-channeling, and high-channeling between rmGAPDH and LDH molecules if we can design experiments that can selectively measure differences in channeling between each of the four subunits in rmGAPDH tetramers. Such experiments will be significantly more demanding than the presented experiments (Supp. Fig. 11), however such efforts can be very significant. The negative cooperativity in NADH binding to different subunits on rmGAPDH was reported almost 60 years ago, yet to this day the physiological significance of such mechanism is not understood 13,47,48 . Differences in NAD(H) binding cooperativity can be observed with GAPDH molecules from different tissues 13,28,46 . Different types of www.nature.com/scientificreports www.nature.com/scientificreports/ binding cooperativity can be also observed between NAD and NADH substrates 13,28,46 . Thus, studies of changes in substrate channeling caused by the changes in NAD(H) binding cooperativity can give crucial insights in the regulation of anabolic and catabolic processes in cells 13,15,17 .

Transient protein-protein interactions and NADH channeling in cells. Presented results (Figs. 4-7)
and the earlier protein-protein interaction studies 20,21 indicate that NADH channeling is taking place within a transient LDH-(GAPDH-NADH) complex (Supp. Fig. 11 and 43 ). Transient interactions could be a physiological necessity for the key metabolic enzymes such as GAPDH which has to rapidly interact with lots of different proteins to accommodate to rapid changes in cell physiology 13,17,27,43 . Transient interactions can be very difficult to describe in studies of protein-protein interaction ( Fig. 7 and 43 ).
The transient GAPDH-LDH complex that forms in experiments with purified enzymes can be much more durable in concentrated protein solution in the cytosol or cell extracts 20,21,51 . The high protein concentrations can produce excluded volume effects that can increase the interaction energy (ΔG binding) by favoring protein-protein interactions that can decrease the number of water molecules trapped in hydration shells 20,21,51 . It is also necessary to notice that in cytosol LDH-GAPDH complex can be a part of much larger glycolytic metabolon 12,15,16 . Such a complex can provide shared molecular scaffolds and molecular crowding effects that can favor collision between the two molecules (on-rates) and decrease the complex breakdown rates (off-rates). Both of those two processes can provide additional stability to LDH-GAPDH complex.
The presented LDH-GAPDH complex indicates that LDH and GAPDH tetramers can simultaneously participate in channeling reaction, in diffusive reaction, and in interaction with other enzymes (supp. Video 1 and Fig. 2B-E). In the presented GAPDH-LDH complex only one of the four subunits in each tetramer is participating in the interaction with its NAD(H) site (Fig. 2C-E). Only NADH binding sites on D subunit on LDH is directly facing the Q subunit on GAPDH (Fig. 2). The NADH binding site on P subunit of GAPDH is open to bind free NAD(H) from the solution, but it is also enclosed in a dominant positive electric field that can channel NAD(H) molecules electrostatically to D subunit on LDH 32 (Fig. 2). A similar situation is observed with the C subunit on LDH (Fig. 2B-E). In other words, free NADH can bind to either Q or P subunits on GAPDH and then get transferred by substrate channeling to D or even C subunits on LDH (Fig. 2B,C). Subunits A and B on LDH and O and R on GAPDH are open for interaction with other enzymes and could not participate in NADH channeling in the presented GAPDH-LDH complex (Supp. Video 1 and Fig. 2B-E). The need to fulfill the multiple functions simultaneously can explain why LDH and GAPDH molecules exist in cells as tetramers. APBS analysis showed that NADH channeling between GAPDH and LDH appears to be an extension of NAD/NADH channeling between the adjacent subunits in GAPDH or LDH tetramers (Figs. 1 and 2).
Substrate channeling is most frequently described as a physiological mechanism that leads to efficient utilization of reactive metabolites 27,33,43 . The concentration of glycolytic enzymes in cells is much higher than the concentration of their substrates 27 , what indicates that in glycolysis the regulation of enzyme activity can be more important than efficient substrate utilization 17,52 . Our results indicate that substrate channeling can regulate metabolism by at least two different mechanisms: (i) allosteric regulation of differences in enzyme-substrate off-rates between different substrates or isozymes (Supp. Baker's yeast cells do not have cytosolic LDH like mammalian cells 53 . Nevertheless, metabolic engineering studies showed that mammalian LDH can readily integrate into glycolic pathways in baker's yeast cells 53 . NADH channeling between byGAPDH and mammalian LDH isozymes can be attributed to highly conserved structures in GAPDH family of enzymes 13 . The conserved structure could be a result of evolutionary pressure to conserve the related interactome 13,17,31 .

Materials and Methods
Apo-enzyme preparations. (NH 4 ) 2 SO 4 suspensions of purified rmGAPDH and byGAPDH have been prepared in our laboratory using established protocols 40,44 or purchased from Sigma Chem. (St. Louis MO).

Molecular dynamics calculations.
All-atom molecular dynamics calculations used the GROMACS 5.1.4 program package 54 (http://www.gromacs.org/) as we have previously described 55,56 . Briefly, NAD(H) molecules from PDB files were processed using ACPYPE10, an interface for Antechamber (part of AmberTools11) that can generate topology types for GAFF force fields 57 . Protein PDB coordinates were processed with pdb2gmx using Amber99SB force field. A cubic solvent box (30 nanometers) was used with TIP3P model for water molecules plus 150 mM NaCl and additional ions that were required for neutralization. The prepared system was minimized using a combination of steepest descent and conjugate gradient algorithms. When the most stable state was achieved the temperature was introduced and the system was equilibrated to 310 K (NVT equilibration, V-rescale). The pressure was equilibrated to 1 atm (NPT equilibration, Parrinello-Rahman). No restraints were used for the protein or the ligand when the system was minimized, but in NPT and NVT equilibration protein and ligand were restrained to prevent the break-up of the complex prior to production runs.
Typical simulations had about 1.5 million atoms, 50 to 150 million steps, with 2 femtosecond time steps. Large simulation boxes (30 nanometers) were used to avoid attractive or repulsive forces created by the periodic boundary conditions. The large boxes can also provide enough space for the two tetramers to dissociate (Supp. Video 1). Different initial simulation set-ups were used to explore optimal design. The simulations have been repeated with the active site loop on LDH in open and closed position, with different initial distances between the interacting proteins, and with different rotations between the LDH and GAPDH tetramers in the interaction plane. Following the simulations, the number of binding interactions was calculated using built-in GROMAC www.nature.com/scientificreports www.nature.com/scientificreports/ functions. The simulations that showed the highest number of binding interactions have been repeated multiple times to assess variability in the number of binding interactions.
Adaptive poisson-boltzmann solver (APBS) calculations. All electric field maps were calculated using Adaptive Poisson-Boltzmann Solver (APBS) approach 58 . PDB formats were converted to PQR format using PDB2PQR and PEOEPB force filed with PROPKA set at pH = 7.2. NAD(H) molecules from PDB files were protonated using GAFF fields. Potential maps were calculated in aqueous 150 mM NaCl solutions using single Debye-Hückel boundary conditions. Fluorescence measurements of NADH binding affinity for rmGAPDH and byGAPDH. Protein fluorescence measurements had excitation 290 nm and emission at 335 nm. NADH fluorescence had excitation at 340 nm and emission at 460 nm 40,44,46,48,50,[59][60][61] . Protein-NADH FRET measurements had excitation 290 nm and emission at 460 nm.
Enzyme buffering measurements. In all enzyme buffering experiments LDH activity was measured by following NADH oxidation in the presence of 630 μM pyruvate which results in a decrease in NADH absorbance at 340 nm. The changes in absorbance were measured with Shimadzu UV-VIS 160 Spectrophotometer or with HP 8452 A Diode Array Spectrophotometer. The molar absorptivity for NADH is 6.22 × 10 3 M −1 cm −1 . All activity measurements were routinely reproducible to with a precession greater than 3%.
The assay mix was prepared in a microcuvette in a total volume of 80 μL and equilibrated to the temperature at 25 °C. The assay buffer was 50 mM Tris/HCl pH = 7.4, 2 mM EDTA-Na, 1 mM DTT, and 1 mg/ml BSA just as in earlier studies 18,19 . First, we filled the cuvette with GAPDH solution in the required concentration and incubated the cell in the holder for several minutes to accommodate to the required temperature and to test the stability of the measured absorbance. Second, NADH was added to the cuvette to confirm that GAPDH solution does not have any intrinsic NADH oxidize activity. Third, pyruvate was added in the concentration of 630 μM next to confirm that there is no measurable LDH activity in GAPDH solution. Finally, the reaction was started by adding 10 to 20 nM of LDH.
In assays with NADH concentration between 40 to 50 μM the initial steady-state rates were measured by following the initial linear decrease in absorbance. In the assays with NADH concentration below 30 μM, the rates were calculated by using exponential equations since NADH concentrations were too low to capture the initial linear decrease in NADH oxidation. byGAPDH solutions were stable even at the highest enzyme concentration tested at 16 mg/ml. rmGAPDH gradually precipitates at the concentration above 8.4 mg/ml, which leads to scattering and detectable increase in absorbance. Thus, the LDH concentrations to maximize the ratio between the changes in absorbance caused by NADH oxidation and by scattering, so that the scattering had to be subtracted as a minor component of the measured changes in NADH absorbance.
The free NADH concentration [NADH]f in the assay mix can be calculated using K D for NADH binding to different subunits on GAPDH (Table 1)

Sedimentation velocity AUC experiments.
An-60-Ti rotor was used with three sample cells with quartz windows and charcoal-epon centerpieces with two sectors 12 mm optical path length. in each cell, 340 μL enzyme samples (6.0 µM) were loaded in one sector and 350 μL buffer in another. Two cells had individual enzymes, and the third cell had the enzyme mixture at the same loading concentration. Three sedimentation velocity experiments were measured in parallel by following absorbance at 280 nm. The first cell had a mixture of phLDH and byGAPDH, both enzymes in the concentration of 6 μM. The second cell had phLDH alone at concentrating of 6 μM, and the third cell had byGAPDH alone in the concentration of 6 μM. The highest protein concentration was determined by the maximal absorbance that can be measured by XL-A instrument (O.D. max = 2.0). The experiments started with the prescans at 3000 × g, and the sum of absorbances of the individual enzymes was compared to the absorbance of the mixture. The sedimentation profiles were measured for 7 hours at 40,000 × g. Single scans with radial increments of 30 μm have been recorded every 240 sec. SEDNTERP was used with protein sequences to calculate molecular mass, partial specific volume, solvent densities, and viscosities 42 . The measured sedimentation profiles were analyzed with SEDFIT program using nonlinear regression with s, and Mr as free fit parameters 42 . www.nature.com/scientificreports www.nature.com/scientificreports/