l-Malate (−2) Protonation State is Required for Efficient Decarboxylation to l-Lactate by the Malolactic Enzyme of Oenococcus oeni

Malolactic fermentation (MLF) is responsible for the decarboxylation of l-malic into lactic acid in most red wines and some white wines. It reduces the acidity of wine, improves flavor complexity and microbiological stability. Despite its industrial interest, the MLF mechanism is not fully understood. The objective of this study was to provide new insights into the role of pH on the binding of malic acid to the malolactic enzyme (MLE) of Oenococcus oeni. To this end, sequence similarity networks and phylogenetic analysis were used to generate an MLE homology model, which was further refined by molecular dynamics simulations. The resulting model, together with quantum polarized ligand docking (QPLD), was used to describe the MLE binding pocket and pose of l-malic acid (MAL) and its l-malate (−1) and (−2) protonation states (MAL− and MAL2−, respectively). MAL2− has the lowest ∆Gbinding, followed by MAL− and MAL, with values of −23.8, −19.6, and −14.6 kJ/mol, respectively, consistent with those obtained by isothermal calorimetry thermodynamic (ITC) assays. Furthermore, molecular dynamics and MM/GBSA results suggest that only MAL2− displays an extended open conformation at the binding pocket, satisfying the geometrical requirements for Mn2+ coordination, a critical component of MLE activity. These results are consistent with the intracellular pH conditions of O. oeni cells—ranging from pH 5.8 to 6.1—where the enzymatic decarboxylation of malate occurs.


Introduction
Most red wines, as well as some white and sparkling wines, are produced by two sequential fermentations: first, yeast alcoholic fermentation (AF) transforms grape must into wine; then, a secondary fermentation, called malolactic fermentation (MLF), is carried out by lactic acid bacteria (LAB). Contrary to AF, a mandatory process of winemaking, MLF is optional and depends on the desired wine style. MLF reduces the acidity of wine and improves flavor complexity and microbiological stability. MLF involves the NAD + -and manganese-dependent decarboxylation of l-malate to l-lactate The malolactic enzyme has been purified from several LAB, e.g., Lactobacillus spp., Lactococcus lactis, and O. oeni [7,17,20,21]. In all cases, the molecular mass of the MLE subunits range from 60 to 70 kDa, and the active form of the protein has been described either as a dimer or a tetramer of identical subunits [22][23][24][25], although all subunits have an independent active site. Interestingly, MLE active sites have binding sites with different structural arrangements of the amino acids Asp, Lys, Ser, and Tyr that, altogether, satisfy the coordination of divalent cation and cofactor positioning [15]. It is noteworthy that all described MLEs catalyze the same reaction.
Among LAB, O. oeni is the best adapted to the harsh conditions of wine. It is capable of carrying out spontaneous fermentation even at pH 3.2, a condition that could be found in some wines [8]. The objective of this study was to provide new insights on the reaction mechanism of the malolactic enzyme (MLE) of O. oeni, responsible for the transformation of L-malate into L-lactate in wine. To this The malolactic enzyme has been purified from several LAB, e.g., Lactobacillus spp., Lactococcus lactis, and O. oeni [7,17,20,21]. In all cases, the molecular mass of the MLE subunits range from 60 to 70 kDa, and the active form of the protein has been described either as a dimer or a tetramer of identical subunits [22][23][24][25], although all subunits have an independent active site. Interestingly, MLE active sites have binding sites with different structural arrangements of the amino acids Asp, Lys, Ser, and Tyr that, altogether, satisfy the coordination of divalent cation and cofactor positioning [15]. It is noteworthy that all described MLEs catalyze the same reaction.
Among LAB, O. oeni is the best adapted to the harsh conditions of wine. It is capable of carrying out spontaneous fermentation even at pH 3.2, a condition that could be found in some wines [8]. The objective of this study was to provide new insights on the reaction mechanism of the malolactic enzyme (MLE) of O. oeni, responsible for the transformation of l-malate into l-lactate in wine. To this end, we first expressed the MLE gene of O. oeni in Escherichia coli BL21. After purification of the protein, we measured the affinity of MLE for malic acid under several pH conditions by isothermal titration calorimetry. Then, we generated a MLE homology model, based on sequence similarity networking and phylogenetic analysis, in order to describe the MLE-malic acid molecular interaction at an atomic level, using molecular dynamics simulations. Finally, we explored the effect of pH on l-malate binding free energies and identified possible residues involved into malic acid binding by means of quantum polarized ligand docking and MM/GBSA calculations. To the best of our knowledge, this is the first study that explores the three-dimensional structure of the malolactic enzyme and its interaction with malic acid at the binding site, the first step of the reaction.

Calorimetric Determination of Malic acid Binding Energies to Malolactic Enzyme
Malolactic fermentation in wine is usually carried out at a pH range between 3.2 and 3.5, allowing a small rise in pH as the malic acid is converted to lactic acid. Figure 2 illustrates the 2D-structure of l-malic acid (MAL) and its l-malate (−1) and (−2) protonation forms (MAL − and MAL 2− , respectively). Isothermal calorimetry thermodynamic (ITC) data for malic acid interaction with MLE shows dissociation constant (K d , Table 1) values in the micromolar range ( Figure 3). MAL 2− has a lower ∆G value than MAL 1− , suggesting this form as the most probable protonation state for malic acid at the binding site of MLE.
end, we first expressed the MLE gene of O. oeni in Escherichia coli BL21. After purification of the protein, we measured the affinity of MLE for malic acid under several pH conditions by isothermal titration calorimetry. Then, we generated a MLE homology model, based on sequence similarity networking and phylogenetic analysis, in order to describe the MLE-malic acid molecular interaction at an atomic level, using molecular dynamics simulations. Finally, we explored the effect of pH on Lmalate binding free energies and identified possible residues involved into malic acid binding by means of quantum polarized ligand docking and MM/GBSA calculations. To the best of our knowledge, this is the first study that explores the three-dimensional structure of the malolactic enzyme and its interaction with malic acid at the binding site, the first step of the reaction.

Calorimetric Determination of Malic acid Binding Energies to Malolactic Enzyme
Malolactic fermentation in wine is usually carried out at a pH range between 3.2 and 3.5, allowing a small rise in pH as the malic acid is converted to lactic acid. Figure 2 illustrates the 2Dstructure of L-malic acid (MAL) and its L-malate (−1) and (−2) protonation forms (MAL − and MAL 2− , respectively). Isothermal calorimetry thermodynamic (ITC) data for malic acid interaction with MLE shows dissociation constant (Kd, Table 1) values in the micromolar range ( Figure 3). MAL 2− has a lower ∆G value than MAL 1− , suggesting this form as the most probable protonation state for malic acid at the binding site of MLE.

Sequence Similarity Networks of Malolactic Enzyme Family
To apprehend the impact of pH on O. oeni MLE activity, we performed an in-silico analysis by comparing its sequence with other MLF-related proteins, including malic enzyme and malate dehydrogenase. For this purpose, we generated a sequence similarity network (SSN), where nodes

Sequence Similarity Networks of Malolactic Enzyme Family
To apprehend the impact of pH on O. oeni MLE activity, we performed an in-silico analysis by comparing its sequence with other MLF-related proteins, including malic enzyme and malate dehydrogenase. For this purpose, we generated a sequence similarity network (SSN), where nodes correspond to homologous proteins to MLE, i.e., those containing at least 70% of sequence identity; and where connections allow to rapidly compute and visualize groups of proteins based on all-against-all sequence comparisons ( Figure 4). The group that contains O. oeni MLE and its closest homologues from NCBI's non-redundant (nr) protein database are referenced as malate dehydrogenases, malic enzymes, and malolactic enzymes. Interestingly, most sequences of this group corresponding to malolactic enzymes and malic enzymes, therefore crystal structures of malic enzymes are the most adequate structural templates to model MLE as there is no structural data available for malolactic enzymes.
against-all sequence comparisons ( Figure 4). The group that contains O. oeni MLE and its closest homologues from NCBI's non-redundant (nr) protein database are referenced as malate dehydrogenases, malic enzymes, and malolactic enzymes. Interestingly, most sequences of this group corresponding to malolactic enzymes and malic enzymes, therefore crystal structures of malic enzymes are the most adequate structural templates to model MLE as there is no structural data available for malolactic enzymes. Clustering by sequence identity is done with CD-HIT program. At values of sequence identity >70%, the nodes should contain sequences that share the same function; however, at lower values of sequence identity, the nodes may be functionally heterogeneous.

Phylogenetics of Malolactic Enzyme Family
Identification of a set of orthologs is a prerequisite for a robust genetic analysis of the evolution of a group of organisms [26]. We carried out multiple sequence alignment using CLUSTAL OMEGA to study the evolutionary relationships between different lactic bacteria in relation to the malolactic enzyme [27]. Most of the sequences homologous to the malolactic enzyme of O. oeni, correspond to proteins whose function has been assigned as malic enzymes by automatic annotation. However, some sequences have been experimentally reported as enzymes with malolactic activity. The latter were labelled with the abbreviation "MLE" below the name of the species, in the phylogenetic tree ( Figure 5).
As illustrated in Figure 5, O. oeni is part of a monophyletic group, together with Streptococcus spp, Lactococcus spp, and Enterococcus spp. However, the evolution of the malic enzyme would be basal in O. oeni with respect to the rest of this clade. It is noteworthy that, within the clade representing the Streptococcus branch, the group of Lactococcus and Enterococcus are represented as sister groups of recent evolution. On the other hand, the branches of Lactobacillus, Pediococcus, and Leuconostoc oeni with at least 70% identity of sequences. The nodes represent proteins and edges indicate similarity in amino acid sequence. Clustering by sequence identity is done with CD-HIT program. At values of sequence identity >70%, the nodes should contain sequences that share the same function; however, at lower values of sequence identity, the nodes may be functionally heterogeneous.

Phylogenetics of Malolactic Enzyme Family
Identification of a set of orthologs is a prerequisite for a robust genetic analysis of the evolution of a group of organisms [26]. We carried out multiple sequence alignment using CLUSTAL OMEGA to study the evolutionary relationships between different lactic bacteria in relation to the malolactic enzyme [27]. Most of the sequences homologous to the malolactic enzyme of O. oeni, correspond to proteins whose function has been assigned as malic enzymes by automatic annotation. However, some sequences have been experimentally reported as enzymes with malolactic activity. The latter were labelled with the abbreviation "MLE" below the name of the species, in the phylogenetic tree ( Figure 5).
As illustrated in Figure 5, O. oeni is part of a monophyletic group, together with Streptococcus spp, Lactococcus spp, and Enterococcus spp. However, the evolution of the malic enzyme would be basal in O. oeni with respect to the rest of this clade. It is noteworthy that, within the clade representing the Streptococcus branch, the group of Lactococcus and Enterococcus are represented as sister groups of recent evolution. On the other hand, the branches of Lactobacillus, Pediococcus, and Leuconostoc constitute a paraphyletic group of basal character with respect to Oenococcus and Streptococcus-Lactococcus-Enterococcus; and they have a previous evolutionary origin.
constitute a paraphyletic group of basal character with respect to Oenococcus and Streptococcus-Lactococcus-Enterococcus; and they have a previous evolutionary origin.

Structural Modeling of the Malolactic Enzyme
The active site of the chain A of the malic enzyme from pigeon liver (PDB entry 1GQ2), the first malic enzyme described [19,28], was selected after SSN analysis as the best structural template for comparative modeling of MLE. The sequence alignment of both structures showed a sequence identity of 35.9% and a coverage of 98% against MLE ( Figure 6). The crystal structure of the A chain contains an oxalate ion in the binding site, and requires Mn 2+ and NADP + as cofactors. Nevertheless, supported by the highly conserved structure of the active site

Structural Modeling of the Malolactic Enzyme
The active site of the chain A of the malic enzyme from pigeon liver (PDB entry 1GQ2), the first malic enzyme described [19,28], was selected after SSN analysis as the best structural template for comparative modeling of MLE. The sequence alignment of both structures showed a sequence identity of 35.9% and a coverage of 98% against MLE ( Figure 6).
constitute a paraphyletic group of basal character with respect to Oenococcus and Streptococcus-Lactococcus-Enterococcus; and they have a previous evolutionary origin.

Structural Modeling of the Malolactic Enzyme
The active site of the chain A of the malic enzyme from pigeon liver (PDB entry 1GQ2), the first malic enzyme described [19,28], was selected after SSN analysis as the best structural template for comparative modeling of MLE. The sequence alignment of both structures showed a sequence identity of 35.9% and a coverage of 98% against MLE ( Figure 6). The crystal structure of the A chain contains an oxalate ion in the binding site, and requires Mn 2+ and NADP + as cofactors. Nevertheless, supported by the highly conserved structure of the active site The crystal structure of the A chain contains an oxalate ion in the binding site, and requires Mn 2+ and NADP + as cofactors. Nevertheless, supported by the highly conserved structure of the active site in both proteins, we confirmed that the putative active binding site could correctly locate malate, after replacing the former cofactors with NAD + and Mn 2+ , and oxalate with malate, using SiteMap of Maestro suite (data not shown). It is worthy to mention that we also employed the malic enzyme from E. coli (PDB entry 6AGS) (Figure 7) for these purposes; though only as secondary scaffold, because no experimental data is available for this protein crystal.
in both proteins, we confirmed that the putative active binding site could correctly locate malate, after replacing the former cofactors with NAD + and Mn 2+ , and oxalate with malate, using SiteMap of Maestro suite (data not shown). It is worthy to mention that we also employed the malic enzyme from E. coli (PDB entry 6AGS) (Figure 7) for these purposes; though only as secondary scaffold, because no experimental data is available for this protein crystal.  Figure 7A shows the MLE homology model we obtained from the abovementioned templates and sequence alignments. This monomeric model was submitted to 200 ns simulation, reaching structural stability after 50 ns, by the structural rearrangement of the carboxyl-term ( Figure S1 in Supplementary materials). Conversely, the pose of NAD + and Mn 2+ reached stability after 20 ns, displaying an RMSD at or below 2 Å throughout simulation. Furthermore, putative MAL bindingsite residues, based on previous reports, namely TYR85, ASP86, LYS156, ASP251, and ASP250 within MLE displayed movement of less than 2 Å ( Figure 7B). Then, MAL was oriented through molecular docking simulations ( Figure 7B).

Molecular Docking of Substrates of Malolactic Enzyme
Additionally, we evaluated the participation of the divalent cation on MLE mechanism by quantum polarized ligand docking (QPLD). Figure 8 illustrates the pose adopted by MAL inside the binding site of MLE. Malic acid interacts with MLE through coordination bonds with Mn 2+ , one LYS protonated residue, and several ASP residues interacting through hydrogen bonds. MAL-MLE interacting residues on this pose correspond with equivalent residues proposed for divalent-cationdependent MAL decarboxylation.  Figure 7A shows the MLE homology model we obtained from the abovementioned templates and sequence alignments. This monomeric model was submitted to 200 ns simulation, reaching structural stability after 50 ns, by the structural rearrangement of the carboxyl-term ( Figure S1 in Supplementary materials). Conversely, the pose of NAD + and Mn 2+ reached stability after 20 ns, displaying an RMSD at or below 2 Å throughout simulation. Furthermore, putative MAL binding-site residues, based on previous reports, namely TYR85, ASP86, LYS156, ASP251, and ASP250 within MLE displayed movement of less than 2 Å ( Figure 7B). Then, MAL was oriented through molecular docking simulations ( Figure 7B).

Molecular Docking of Substrates of Malolactic Enzyme
Additionally, we evaluated the participation of the divalent cation on MLE mechanism by quantum polarized ligand docking (QPLD). Figure 8 illustrates the pose adopted by MAL inside the binding site of MLE. Malic acid interacts with MLE through coordination bonds with Mn 2+ , one LYS protonated residue, and several ASP residues interacting through hydrogen bonds. MAL-MLE interacting residues on this pose correspond with equivalent residues proposed for divalent-cation-dependent MAL decarboxylation. We also calculated the most probable protonation state of MAL using the Epik module of the Schrodinger Suite. Results confirmed that MAL 2− is the most probable protonation state and thus interacting residues could be oriented differently to MALand MAL.
Quantum polarized ligand docking (QPLD) was then employed to explore MAL 2− pose and binding energy (∆Gbinding). Results showed that all malic acid protonation forms lie in the same binding cavity, sharing the same set of binding amino acid residues ( Figure 8). The latter interact mainly by hydrogen bonding and hydrophobic interactions ( Figure S2 in Supplementary materials). MAL 2− has the lowest ∆Gbinding, followed by MAL − and MAL, with values of −23.8, −19.6, and −14.6 kJ/mol, respectively (Table 2). Interestingly, MAL 2− displayed an extended conformation, when compared with the other protonation states. This open conformation better satisfies the geometrical requirements of Mn 2+ coordination geometry and the mechanism described for malic enzymes [29]. We also calculated the most probable protonation state of MAL using the Epik module of the Schrodinger Suite. Results confirmed that MAL 2− is the most probable protonation state and thus interacting residues could be oriented differently to MALand MAL.
Quantum polarized ligand docking (QPLD) was then employed to explore MAL 2− pose and binding energy (∆G binding ). Results showed that all malic acid protonation forms lie in the same binding cavity, sharing the same set of binding amino acid residues ( Figure 8). The latter interact mainly by hydrogen bonding and hydrophobic interactions ( Figure S2 in Supplementary materials). MAL 2− has the lowest ∆G binding , followed by MAL − and MAL, with values of −23.8, −19.6, and −14.6 kJ/mol, respectively (Table 2). Interestingly, MAL 2− displayed an extended conformation, when compared with the other protonation states. This open conformation better satisfies the geometrical requirements of Mn 2+ coordination geometry and the mechanism described for malic enzymes [29].
Geometrical stability of MAL inside the binding site was assessed after 200 ns molecular dynamics simulations of the MLE/MAL/NAD + /Mn 2+ system. Of note, MAL does not remain on the site and exits the pocket at 25 ns. On the contrary, MAL − and MAL 2− remain into the binding pocket throughout the whole simulation. Moreover, the average binding energy of the molecular interactions through the MM/GBSA rescoring method was calculated as it is relatively more accurate compared to single-structure theoretical determinations. MM-GBSA binding energies for MAL − and MAL 2− showed binding affinity differences consistent with values from the ITC and QLPD experiments (Tables 1 and 2, respectively). Furthermore, energy decomposition of MAL − and MAL 2− interactions within the MLE binding pocket allowed to identify that binding is mainly driven by the negative charge interactions of the MAL carboxyl group with the positive charge of the side chain-N of LYS156, (Figure 7); whereas MAL 1− interacts with ASP86 and ASN396 mainly by H-bonding; and MAL 2− with TYR85, ASP86, ASP228, ASN396, and ASN440 mainly through water bridges.     Geometrical stability of MAL inside the binding site was assessed after 200 ns molecular dynamics simulations of the MLE/MAL/NAD + /Mn 2+ system. Of note, MAL does not remain on the site and exits the pocket at 25 ns. On the contrary, MAL − and MAL 2− remain into the binding pocket throughout the whole simulation. Moreover, the average binding energy of the molecular interactions through the MM/GBSA rescoring method was calculated as it is relatively more accurate compared to single-structure theoretical determinations. MM-GBSA binding energies for MAL − and MAL 2− showed binding affinity differences consistent with values from the ITC and QLPD  Geometrical stability of MAL inside the binding site was assessed after 200 ns molecular dynamics simulations of the MLE/MAL/NAD + /Mn 2+ system. Of note, MAL does not remain on the site and exits the pocket at 25 ns. On the contrary, MAL − and MAL 2− remain into the binding pocket throughout the whole simulation. Moreover, the average binding energy of the molecular interactions through the MM/GBSA rescoring method was calculated as it is relatively more accurate compared to single-structure theoretical determinations. MM-GBSA binding energies for MAL − and MAL 2− showed binding affinity differences consistent with values from the ITC and QLPD  Geometrical stability of MAL inside the binding site was assessed after 200 ns molecular dynamics simulations of the MLE/MAL/NAD + /Mn 2+ system. Of note, MAL does not remain on the site and exits the pocket at 25 ns. On the contrary, MAL − and MAL 2− remain into the binding pocket throughout the whole simulation. Moreover, the average binding energy of the molecular interactions through the MM/GBSA rescoring method was calculated as it is relatively more accurate compared to single-structure theoretical determinations. MM-GBSA binding energies for MAL − and MAL 2− showed binding affinity differences consistent with values from the ITC and QLPD

Discussion
In silico analysis, MLE homology model together with QPLD and ITC experiments were carried out in the present study to determine the protonation state of l-malate required for its efficient decarboxylation to l-lactate by the malolactic enzyme of Oenococcus oeni.
The phylogenetic analysis showed that the evolution of the O. oeni malolactic enzyme is halfway between malic-malolactic enzymes of the genera Lactobacillus-Pediococcus-Leuconostoc and Streptococcus-Lactococcus-Enterococcus. The grouping of lactic acid bacteria in two clusters was in line with Makarova and Koonin (2007) [30], which employed the sequence of ribosomal proteins and RNA polymerase subunits for their phylogenetic analysis. Our results indicated that O. oeni MLE is part of a monophyletic group, together with the branch Streptococcus-Lactococcus-Enterococcus, whereas O. Oeni MLE constitutes a paraphyletic group of the Lactobacillus-Pediococcus branch. These results point that the genera Streptococcus, Lactococcus, Enterococcus, and Oenococcus descended from a common evolutionary ancestor, whose malolactic enzyme could share similar structural characteristics.
Malolactic fermentation usually occurs at pH levels between the range of 3.2 and 3.5, allowing for a rise in pH as the malic acid is converted to lactic acid. At this pH, the MAL − protonated l-malic acid form prevails (Figure 2). On the other hand, several LAB, including O. oeni have an intracellular pH ≈ 5.0 [31,32], a condition where MAL 2− is the predominant protonated form. Several studies have reported that the malolactic fermentation reaction occurs in the Oenococcus oeni intracellular space, where pH is between 5.8 and 6.1; while in the extracellular medium, i.e., the wine conditions, the pH is within a range of 3 to 4 [33,34]. Additionally, Schümann et al. (2013) reported that O. oeni MLE has an optimum activity at pH 6.0 and 45 • C [14]. Accordingly, we evaluated the effect of pH on l-malic acid interaction with the MLE active site. To this end, we measured the enthalpy of reaction of O. oeni MLE with both cofactors, titrated with MAL at pH 4.5 and 5.3, using ITC. Our results showed a higher binding affinity for MAL 2− than for MAL − , in agreement with Schümann et al. (2013) [14]. Under the O. oeni intracellular conditions, the presence of some residues, in the binding pocket or in its vicinity, could accept the proton from the MAL − form, predominant in solution, to lead the most stable MAL 2− , such as Asp86 and Glu227 (see Figure 8A).
Although ITC experiments showed that pH significantly influences malic acid binding to O. oeni MLE, this method did not allow to extract structural information of the binding sites or enzymatic mechanisms, at atomic level. However, non-integer stoichiometric values indicated the formation or aggregation of dimers of higher quaternary structures, which was also observed by Dynamic Light Scattering measurements (data not shown). These results agree with the work of Schümann et al. (2012), where MLE of O. oeni was presented as a dimeric macromolecule, with each subunit having a functional binding site [35].
To give further structural insights and to understand the calorimetric results, we relied on the use of molecular docking and molecular dynamics methods. To this end, O. oeni MLE homology model was carried out, because the three-dimensional structure was not available. It is worth noting that SSN and phylogenetic analysis showed a close relation between malolactic and malic enzymes (Figure 4), although only four crystal structures are available as possible structural templates. Among these structures, malic enzyme from pigeon liver (PDB entry 1GQ2) and E. coli Malic enzyme (PDB entry 6AGS) were used as model templates, using the alignment shown in Figure 6. Both structures were identified as the closest related sequences according to the SSN; nevertheless, experimental data regarding 6AGS crystal is scarce, and was only used when no structural data from 1GQ2 were available. It is worth noting that NAD + , Mn 2+ were incorporated as model constraints, using the pose of cofactors found on the 1GQ2 crystal. The MLE model was submitted to 200 ns molecular dynamics simulations to further explore its dynamic and stability and to identify relevant residues for malic acid binding. The trajectory analysis showed that, overall, residues near cofactors have slow mobility and form cavities that are suitable to bind l-malic acid.
Interestingly, SSN and phylogenetic analysis suggest a closed relation between malate dehydrogenases, malic enzymes, and malolactic enzymes; however, their binding sites are not completely conserved and thus, substrate pose, and binding residues could be oriented differently. To correctly orientate malic acid and to calculate theoretical binding energies considering Mn 2+ , we opted for the quantum polarized ligand docking (QPLD) method as it allows to properly calculate binding energies of metal-containing systems as it considers metal coordination, electronic polarization effect, among other missed terms in molecular-mechanics force fields. As can be seen in Table 2, QPLD results shown that MAL − and MAL 2− have docking scores (∆G binding ) that correlate with ITC measurement results. Although both MAL protonation states bind into the same cavity and share a set of amino acids, MAL 2− adopts an extended conformation, supported by hydrogen bonding and coordination geometry with Mn 2+ (Figure 7). Further, we explore the binding site dynamics of the MAL, MAL − , and MAL 2− containing systems through molecular dynamics simulation. Of note, MAL 2− kept its orientation, as determined by the QLPD method, after 200 ns of molecular dynamics simulation, while MAL − and MAL have more mobility inside the binding site. Furthermore, molecular dynamics provides a conformational ensemble that allows to calculate the average binding energy of the molecular interactions through the MM/GBSA rescoring method that is more accurate compared to single-structure theoretical determinations as it includes solvent effects. According to MM/GBSA results and ITC experiments, MAL 2− have the lowest binding energy (∆G binding ) and the major energetic contribution is the stabilization of the two carboxylic group charges that interacts with Mn 2+ . Regarding ITC correlation with our molecular dynamics results, it should be noted that docking and MMGBSA calculations represent one of the steps of the reaction coordinates, that is pre-transition states without consider the diffusion pathways into the binding sites and omitting desolvation and other effects that directly impact the entropy variations observed by ITC measurements. This MAL 2− pose is in accordance with the mechanism described by Schümann et al. (2013), where Mn 2+ act as an activator of the enzymatic catalysis and coordinate chemical reaction, while NAD + act as oxidizing agent for oxidation of l-malate.

Analysis of Sequences and Construction of Phylogenetic Tree
The amino acid sequence of malolactic enzyme (MLE) from Oenococcus Oeni was retrieved from NCBI (accession number WP_002823502.1). Homologous amino acid sequences were found using the online available version of Basic Local Alignment Tools (BLAST™) [36]. 1, and ARC49389.1) of lactic acid bacteria were selected based on the e-value, the query coverage and its sequence identity with MLE, and were aligned using CLUSTAL OMEGA program (EMBL) [27]. The phylogenetic tree was built up using Interactive Tree Of Life [37].
To perform the Sequence similarity Network, the MLE sequence from O. oeni (Uniprot ID: Q48796) was aligned to the closest sequences (>70% id) using all-by-all BLAST within InterProScan database performed by the web service EFI-Enzyme Similarity Tool [38,39]. Each sequence was labeled by its primary biological function and structural data availability as provided by UNIPROT. Even though all the function entries were cured manually, the lack of consistency of UNIPROT terminology could lead to ambiguous descriptions. Finally, 309,755 sequences were admitted to the SSN building.

Template Selection
The amino acid sequences of the malolactic enzyme from O. oeni strain DSM 20255 were retrieved from the NCBI (accession number ACX50963). The template was selected based on the e-value of the BLAST search, query coverage, and its sequence identity with MLE. Based on these criteria, malic enzyme from pigeon liver (PDB entry 1GQ2) was selected as a template to model MLE.

Modeling of Malolactic Enzyme
A comparative model for the malolactic enzyme was constructed using Prime from Schrodinger Suite 2019 using the PDB 1GQ2 and 6AGS as a template. Both enzyme cofactors, NAD + and Mn +2 , were incorporated into the resulting models keeping the atomic coordinates from reference structure 1GQ2. The resulting model of the malolactic enzyme, including NAD + cofactor and Mn +2 ion, was inserted in a water box and further neutralized with counter ions. Then MLE/NAD + /Mn +2 complex was subjected to cycles of energy minimize as described elsewhere and equilibration for 200 nanoseconds under NPT conditions.

Ligand Preparation
l-malic acid three-dimensional structure was obtained from the PubChem database (Pubmed CID 222656) and prepared in Maestro (Schrödinger, LLC, New York, NY, USA) using the OPLS_2005 force field with default setting of the LigPrep package from Schrödinger. All molecules were visualized and pKa values were calculated using Epik at desire pH [40].

Quantum Polarized Ligand Docking (QPLD)
l-malic acid and its other protonation states were docked with improved docking program of quantum polarized ligand docking (QPLD) of the Schrodinger Suite 2019 [41]. The best poses obtained by flexible ligand docking using Glide [42]. Then QM calculations were done using Jaguar to calculate the partial charges were replaced on the ligand in the field of receptor for each ligand complex [43]. Single point electrostatic calculations were carried out with the 6-31G*/LACVP* base set and B3LYP density functional theory, using the "Ultrafine" SCF accuracy level (iacc = 1, iac-scf = 2) for the QM region. Finally, ligand was redocked with updated atomic charges with the help of Glide XP and QPLD of the Schrodinger Suite 2019.

Molecular Dynamics Simulation (MD)
MD simulations of malolactic enzyme and ligand complexes were carried out using Desmond and OPLS 2005 force field [44]. The protein ligand complexes were solvated with TIP3 water molecules. Sodium counter ions were added to balance the system net charge. The systems were submitted to the default Desmond protocol, which contains a series of restrained minimizations and MD simulations. The minimized system was relaxed under NPT ensemble for 50 ns equilibration simulation period, and 150 ns production simulations were carried out. Long range electrostatic interactions were computed by particle-mesh Ewald method and van der waals (VDW) cut-off was set to 9 Å.
Escherichia coli BL21 strain and plasmid pET28a were obtained from Novagen (Buenos Aires, Argentina). Transformants were grown at 37 • C in LB medium, with the addition of 50 µg/mL kanamycin. Agar plates were made of LB media, including 15 g/l agar.
All PCRs to amplify DNA fragments suitable for Gibson assembly were carried out in 35 PCR cycles, using Phusion High-Fidelity DNA polymerase (Thermo Scientific, Waltham, MA, USA), following the manufacturer's instructions. Gibson assembly was performed as previously described [48] with pairs of primers for each fragment to be assembled containing segments of about ∼40 bp homologous to the adjacent fragment to be linked. All PCR products were treated with the DpnI enzyme to eliminate original vector residues and purified by gel extraction using the Qiaquick Gel Extraction kit from Qiagen, according to the manufacturer's instructions. The purified genes fragments and vectors were mixed based on their molar ratios in a final volume of 5 µL, containing 100 ng of total DNA. This DNA mix was added to 15 µL of 1.33X master mix (5X isothermal mix buffer, T5 exonuclease 1 U/µL, Phusion DNA polymerase 2 U/µL, Taq DNA ligase 40 U/µL and Milli-Q purified water), and the reaction mixture was incubated at 50 • C for 1 h. Finally, 10 µL reaction mix were used directly to transform chemically competent E. coli BL21 (DE3).
The vector construct, designated pET21a-MLE, was verified by sequencing (Macrogen Inc., Seoul, Korea). The resulting map is shown in Figure S3  The resulting biomass was recovered from the fermentation broth by centrifugation (4000x g, 10 min, 4 • C) and the supernatant was discarded. Approximately 9 g of biomass were recovered from 1 L of fermentation broth. Subsequent cell disintegration was carried out in lysis buffer (Tris 20 mM pH 6.0, with 500 mM NaCl, 30 mM imidazole, and protease inhibitor cocktail complete™), at a concentration of 1 g of biomass in 10 mL of lysis buffer. The mix was distributed in 1.5 mL Eppendorf tubes with 250 µL of glass beads (Sigma-Aldrich ® ), and cell disruption was performed by agitation, three consecutive cycles of 30 s.

Calorimetric Characterization
Enthalpy changes associated with MLE-substrate interactions were measured using a Nano ITC instrument (TA Instruments Ltd., Crawley, West Sussex, U.K.), at 25 • C. An amount of 170 µL of MLE solution (30 µM, HEPES buffer at desired pH) were placed in the sample cell of the calorimeter and buffered substrate solution (100 µM, HEPES buffer at desired pH) was loaded into the injection syringe. The substrates were titrated into the sample cell as a sequence of 20 injections of 2.5 µL aliquots. The time delay (to allow equilibration) between successive injections was 3 min. The contents of the sample cell were stirred throughout the experiment at 200 rpm to ensure thorough mixing. Raw data were obtained as a plot of heat (µJ) against injection number and featured a series of peaks for each injection. These raw data peaks were transformed using the instrument software Nano Analyze (version 3.11.0, TA Instruments, New Castle, DE, USA) to obtain a plot of observed enthalpy change per mole of injectant against molar ratio and were corrected by subtracting the mixing enthalpies of the substrate solutions into protein-free solution.

Conclusions
In conclusion, in this work, we constructed a comparative model for MLE using the 3D structures of the malic enzyme from pigeon liver (PDB entry 1GQ2) and malic enzyme from E. coli (PDB entry 6AGS). Malic acid interactions within MLE binding pocket are mainly driven by hydrogen bonding and coordination with Mn 2+ , both dependent on the protonation state of the substrate.
Our experimental and theoretical studies demonstrated that MAL 2− stabilizes the pose that fulfills the geometrical requirements to favor the malic acid decarboxylation catalyzed by MLE. Further theoretical and experimental studies are currently underway to provide more detailed information about the contribution of each residue on the MLE proposed mechanism.
Supplementary Materials: The following are available online. Figure S1: