Substrate Specificity of the Highly Thermostable Esterase EstDZ3

Esterases are among the most studied enzymes, and their applications expand into several branches of industrial biotechnology. Yet, despite the fact that information on their substrate specificity is crucial for selecting or designing the best fitted biocatalyst for the desired application, it cannot be predicted from their amino acid sequence. In this work, we studied the substrate scope of the newly discovered hydrolytic extremozyme, EstDZ3, against a library of esters with variable carbon chain lengths in an effort to understand the crucial amino acids for the substrate selectivity of this enzyme. EstDZ3 appears to be active against a wide range of esters with high selectivity towards medium‐ to long‐carbon chain vinyl esters. In‐silico studies of its 3D structure revealed that the selectivity might arise from the mainly hydrophobic nature of the active site's environment.


Introduction
Lipolytic enzymes have received a great deal of attention as they are ubiquitous in nature and are easily observed in all kingdoms of life. Carboxylic ester hydrolases (EC 3.1.1.x) represent a broad class of lipolytic enzymes catalysing the hydrolysis of ester bonds over a wide range of substrates. This class of enzymes can be further subcategorized into several groups based on their substrate specificity. Among them, two subclasses, carboxylesterases (EC 3.1.1.1) and triacylglycerol hydrolases (EC 3.1.1.3, also known as lipases), have been excessively studied. [1] There have been many attempts from the scientific community to establish an adequate distinction between these classes through the establishment of various criteria. Lipases and esterases act upon water-insoluble and -soluble substrates, respectively. Moreover, a distinctive characteristic of lipases is that they can be interfacially activated. [2] This characteristic is attributed to a "lid" structure, composed from an amphiphilic α-helix that covers the hydrophobic substrate pocket. [3] However, detailed biochemical and structural studies have proven these criteria to be not exact determinants. For example, some lipases exhibit high affinity for water-insoluble esters without possessing a lid domain. [4] Lipolytic enzymes are of great biological importance due to the fact that they contribute to a variety of biological processes [5] Indeed, their functions are critical for human physiology owning to the fact that they can provide carbon sources for energy, through triacylglycerol catabolism, while the resulting products can also act as precursors or mediators for various biosynthetic and cellular signalling processes. [6,7] Some of their most industrially favourable characteristics include stability in organic solvents, their ability to catalyse ester synthesis in unconventional media, as well as their capacity to catalyse reactions with high chemo-, regio-, and enantioselectivity without the need of cofactors. All these features render them greatly important catalysts for industrial biotransformations. [8,9] Some notable industrial applications include the production of dairy products, the degradation of plastics, the synthesis of fine chemicals and their use as diagnostic tools. [10][11][12][13] More recently, the focus is gradually shifting towards the development of novel immobilization carriers to improve the stability of such biocatalysts. An example of such an application recently published [14] has clearly showcased the advantages of an ecologically friendly process based on the immobilization of a Thermomyces lanuginοsus lipase on hybrid magnetic ZnOFe nanoparticles, which greatly enhanced the stability of the enzyme leading to high transformation yields.
The vast majority of lipolytic enzymes incorporated in industrial settings are of microbial origin. For this reason, there is a growing interest in discovering new enzymes with enhanced thermal stability, which can withstand high-temperature industrial processes. Such enzymes are frequently encountered in organisms residing in extreme environments. [15] One powerful strategy to identify such biocatalysts is metagenomic analysis, either by bioinformatic or functional means. [16][17][18] Metagenomics offer the great advantage that they can bypass the limitations imposed by traditional microbial culturing techniques. [19] We have recently discovered the new esterolytic enzyme EstDZ3, which was found to exhibit remarkable thermostability. [20] EstDZ3 was isolated from a bacterium of the genus Dictyoglomus, identified in a hot spring in China. [20] Interestingly, this new enzyme was found to have low amino acid sequence similarity to any previously identified enzyme. Molecular modelling of its three-dimensional structure has provided an indication of the existence of a "subdomain insertion", which is similar to that of the cinnamoyl esterase Lj0536 from Lactobacillus johnsonii. [20,21] According to our initial biochemical studies, EstDZ3 seems to function more like a carboxylesterase rather than a lipase, as it demonstrated a clear preference towards p-nitrophenyl (pNP) esters with short to medium acyl-chain length [20] More importantly, EstDZ3 was stable when exposed to high concentrations of organic solvents and exhibited remarkable thermostability, characteristics which render it a promising new catalyst for industrial biotransformations.
In this work, we have investigated in more detail the hydrolytic activity and specificity of EstDZ3 against a variety of aliphatic and aromatic esters of synthetic or natural origin. Furthermore, we have performed computational analyses to provide rationalization for the identified substrate specificity. Our studies contribute important understanding of the biochemical properties and the biocatalytic potential of this new enzyme.

Results and Discussion
The esterolytic activity of EstDZ3 We investigated the hydrolytic activity of EstDZ3 against various synthetic and natural aromatic esters ( Figure S1 in the Supporting Information). As seen in Figure 1, EstDZ3 can hydrolyse a broad range of aliphatic esters, as well as aromatic esters, albeit with a lower efficiency. Among the aromatic esters, EstDZ3 exhibited good activity against benzyl acetate and cinnamyl acetate, while low activity was observed against bulky natural substrates, such as oleuropein and rosmarinic acid. On the other hand, the enzyme exhibited up to 20-fold higher hydrolytic activity towards vinyl esters ( Figure 1). This observation may be attributed to the fact that the electron-withdrawing groups of vinyl esters are able to shift the reaction equilibrium towards hydrolysis through the instant tautomerization of vinyl alcohol to the highly volatile acetaldehyde. [22] The effect of the acyl chain length (from C 2 to C 12 ) on the hydrolytic activity of EstDZ3 was investigated using aliphatic vinyl esters. EstDZ3 showed a preference towards medium to long aliphatic esters (C 10 À C 12 , Figure 1b highest hydrolytic activity towards vinyl decanoate (C 10 ), while lower levels of activity were observed for substrates with shorter or longer carbon chains.
A more extensive and detailed study of EstDZ3's selectivity was conducted by determining its apparent kinetic parameters towards the aforementioned aliphatic esters through the Michaelis-Menten kinetic model (Table 1). EstDZ3 displayed the highest catalytic efficiency for vinyl decanoate with a k cat /K m value of 3511 s À 1 mM À 1 , while its ability to hydrolyse esters with a carbon chain length longer or shorter than C 10 was significantly lower. The lowest esterolytic activity was observed against vinyl acetate, with a k cat /K m value of 95 s À 1 mM À 1 . Based on these findings, it is evident that EstDZ3 is not a typical esterase as it accepts larger substrates than other esterases.
It is interesting to note that these results are not in accordance with our initial substrate specificity analysis using pNP esters. [20] In that work, EstDZ3 displayed a high preference towards pNP esters with short to medium carbon chain length (C 2 À C 8 ), as expected for a typical esterase, with the highest affinity observed for p-nitrophenyl butyrate. The different specificity observed in the present work indicates that the overall size of the substrate might affect the affinity of the enzyme for each substrate. pNP esters are bulkier substrates compared to vinyl esters. Moreover, they differ in the alcohol moiety, which in turn might also have a potential effect on the overall specificity of EstDZ3.

Docking esters into the homology model of EstDZ3
We continued by employing computational methods in an attempt to gain a deeper structural insight and rationalize the substrate specificity of EstDZ3. Currently there is no available three-dimensional structure for EstDZ3. The published EstDZ3 homology model, as predicted by iTasser [23] revealed that its 3D structure is characterized by the typical α/β hydrolase fold, while the catalytic triad is comprised by the residues S114, D202 and H233 [20] In order to obtain a refined 3D structure of EstDZ3 for our in silico studies, we utilized the algorithm provided by AlphaFold. [24] The 3D structure of the protein can be found in Figure S2.
Molecular docking was performed to elucidate the binding mode of vinyl esters, for which kinetic data were obtained, to EstDZ3 and to provide an insight into the defining elements of the enzyme specificity. Table 2 summarizes the docking results (productive binding poses) of EstDZ3 against various vinyl esters with variable acyl chain sizes. As observed in the table, all substrates were docked successfully into the binding pocket of EstDZ3 with binding affinities ranging from À 3.03 to À .31 kcal mol À 1 . Vinyl decanoate exhibited the lowest binding energy as well as the lowest K d , which indicates that this substrate binds tightly into the active site. Regarding vinyl acetate, it is apparent that the acetate group cannot fully exploit the binding pocket for interactions (Table 2, Figure S3e), while vinyl decanoate (C 10 ) seems to be the optimal substrate size ( Figure S3b). Vinyl laurate (C 12 ) is positioned in a way which yields a bent of the alkyl group towards the solvent (Figure S3a), indicating that the size of the substrate is probably larger than what can be accommodated in the pocket. These results are in agreement with experimental K m values. It seems that active site of EstDZ3 cannot fully accommodate substrates with C 12 atoms or larger, hence the observed sharp decrease in enzyme activity for longer substrates. In addition, the binding  mode for almost all substrates is similar. More specifically, the vinyl moiety of the substrate positions itself closer to the active site residues, S114 and H233, while the acyl chain expands and adopts a proper orientation into the rest of the enzyme's binding pocket that seems to be favoured by hydrophobic interactions. These docking results are in accordance with the experimental results described in the previous paragraph, as EstDZ3 exhibited the highest affinity for vinyl decanoate and the lower affinity for vinyl acetate. It is noteworthy that the environment of the enzyme's active site is composed mainly by hydrophobic amino acids such as valine, glycine, and leucine ( Figure S4). The important role of the hydrophobic interactions is also highlighted in Table 2, since most of the interactions recorded are between hydrophobic residues, apart from the polar interactions with the catalytic residues and the oxyanion hole. This might serve as a possible explanation on why the selectivity of the enzyme is towards aliphatic substrates, that is, more hydrophobic, instead of aromatic ones, which indicates a better stabilization inside the active site. Furthermore, as already mentioned above, the closest homolog to EstDZ3 that has been structurally and biochemically characterized, is the cinnamoyl esterase Lj0536. This protein possesses an "inserted sub-domain" in close proximity to the active site, [21] which appears to play an important role for the ability of the enzyme to discriminate among different substrates. The sub-domain of Lj0536 is mainly comprised of polar residues ( Figure S4b) as opposed to that of EstDZ3, which is mainly composed of hydrophobic ones ( Figure S4a). This difference might also explain why Lj0536 displays higher preference towards aromatic esters [21] as opposed to the aliphatic ones showcased for EstDZ3 in this study.

Molecular dynamics simulation studies
Molecular dynamics simulations were performed on substrates with different alkyl-chain lengths, to validate the stability of binding modes of the vinyl esters. For this reason, we focused on vinyl decanoate (C 10 ) and vinyl butyrate (C 4 ), which differ significantly in terms of binding energy and acyl-chain length. In order to investigate the ability of the candidate substrate to maintain a proper catalytic orientation, we measured the distances of the critically relevant atoms between the carbonyl carbon atom (Ca) of the substrates and the Oγ of the catalytic serine as well as the distances between the carbonyl oxygen and the amide group of the oxyanion hole residues as a function of time during the simulation (Figures S5 and S6). In Figure 2, the distances between the substrates' Ca and the Oγ atom in catalytic serine are represented as a function of simulation time. In the case of vinyl butyrate (Figure 2a), the distances drastically increase after 10 ns, ranging from~10 tõ 67 Å, indicating that the substrate is no longer accommodated in a catalytic orientation inside the active site, and instead dissociates in the solvent. Indeed, as seen in Figure 3a, vinyl butyrate occupies a different position by the end of the simulation way far from the protein's active site. This orientation seems to be favoured by the hydrophilic interactions between the substrate's carbonyl oxygen and water molecules of the simulation's solvent (molecules not shown). Therefore, vinyl butyrate fails to maintain a stable conformation close to the catalytic residues due to the lack of hydrophobic interactions between the substrate's short acyl chain and the enzyme's active site. These results indicate a low affinity of the enzyme for vinyl butyrate. Regarding vinyl decanoate, the results from the molecular dynamics analysis have shown a different outcome. It is observed that, at the beginning of the simulation, the carbonyl carbon maintains a distance smaller than 4 Å from serine's Ογ atom. Throughout the course of the simulated time, fluctuations arise on the substrate's conformation, affecting its orientation in the active site (Figure 2b). The substrate seems to be moving further from the catalytic serine, which is evident from the increase of its distance, but also seems to maintain its catalytic orientation (< 4 Å) in the majority of the spectra's simulated time (Table S2). This is a great indication that in the timeframe of the analysis the substrate preserved a catalytic conformation. Additionally, the distances between the carbonyl oxygen and the residues of oxyanion hole (Phe38 and Met115) represent a similar behaviour ( Figure S6). In Figure 3b, it is evident that the orientation of the substrate mainly in the end of the simulation, involves the acyl-chain moiety and the vinyl moiety as well. This movement of the acyl-chain is likely because of the core hydrophobic pocket near the active site, which seems to be favoured by the hydrophobic interactions. Nevertheless, in Table S2 it is still evident that for the majority of time, the substrate remained in a catalytically active conformation based on the critical distances between the atoms of the complex.
The aforementioned observations from the molecular dynamics analysis are in accordance with our experimental results and also further confirm our docking analysis. Vinyl decanoate is able to maintain a more stable and catalytically active orientation inside the active site, which indicates that the enzyme displays a higher affinity for this substrate and, thus, can exhibit higher catalytic activity against it compared to vinyl butyrate.

Conclusions
In this work, we have showcased the ability of EstDZ3 to hydrolyse an extended spectrum of substrates. EstDZ3 was able to hydrolyse vinyl esters with great catalytic efficiency, with maximum efficiency against vinyl decanoate (C 10 ). EstDZ3 exhibited lower activity towards bulky and aromatic substrates, a characteristic that may be attributed to the fact that its active site is mainly hydrophobic. Docking analysis revealed that the esterolytic specificity is a factor of the chemical environment of the enzyme's active site, as it was able to accommodate substrates with carbon chain lengths up to C 10 with great ease. At the same time, hydrophobic interactions were defining for the stability of these substrates in the binding pocket. Further analysis from molecular dynamics simulations validated our docking results, which in turn were in accordance with our experimental observations. In conclusion, EstDZ3 specificity is affected not only by the substrates' fatty acid chain length, but also it is strongly dependent on the nature of the alcohol part of the ester substrates.

Experimental Section
All chemical reagents were purchased form Sigma-Aldrich or Fluka. All organic solvents were HPLC grade.
Enzyme expression and purification: EstDZ3 recombinant production and purification was carried out as described previously. [20] Briefly, BL21(DE3) cells carrying the pLATE52-EstDZ3 plasmid were grown at 37°C until the culture reached an optical density at 600 nm (OD 600 ) of about 0.5. At that point, the expression of estDZ3 was induced by the addition of isopropyl-β-d-thiogalactoside (IPTG) followed by overnight incubation at 25°C. The cells were harvested and EstDZ3 was recovered using the Qiagen IMAC purification kit according to the manufacturers' instructions with the following modifications. The cell pellet was washed and resuspended in equilibration buffer supplemented with 1 % TritonX-100, and lysed by brief sonication. The cell extract was clarified by centrifugation and the supernatant was incubated for 30 min at 80°C in order to denature Escherichia coli-soluble proteins. After heat-treatment, the precipitated material was removed by centrifugation. The supernatant was collected and incubated with Ni-NTA agarose beads before it was loaded onto a polypropylene column (ThermoScientific). The flow-through was discarded, and the column was washed with NPI20 wash buffer containing 1 % Triton X-100. Next, Triton X-100 was washed away by passing standard NPI20 wash buffer through the column. EstDZ3 was eluted using NPI200 elution buffer. Imidazole was subsequently removed from the protein preparation using a Sephadex G-25 M PD10 column (GE Healthcare). EstDZ3 was retrieved in apparent purity according to SDS-PAGE analysis visualization. The protein preparation was lyophilized and stored at À 20°C.
Enzyme activity assays: A colorimetric method was used for the determination of the substrate profile of EstDZ3, [25,26] with slight modifications. The purified and lyophilized EstDZ3 was dialysed in BES [N,N-bis(2-hydroxyethyl)-2-aminoethanesulfonic acid] buffer 2.5 mM, pH 7.2 at 1 mg mL À 1 . All substrates were diluted in acetonitrile at 30 mM. A standard reaction mixture contained 0.23 mM of the pH indicator, p-nitrophenol (pNP), 7.1 % acetonitrile, 1 mM substrate and 30 μg mL À 1 enzyme EstDZ3. The enzymatic reaction was performed in a total volume of 1050 μL in a cuvette, in an Agilent Technologies Cary 60 UV-Vis spectrophotometer, equipped with a Peltier system for temperature control. The decrease in absorbance was monitored at 405 nm, at 50°C, for 10 min, at 11 s intervals. The extinction coefficient of pNP at the above conditions was found to be ɛ = 10 100 M À 1 cm À 1 .
To determine the kinetic parameters K M and V max , an enzyme concentration of 5 μg mL À 1 was used, and the substrate concentration ranged from 0.02 to 2 mM. The EnzFitter Biosoft (UK) software was used for data analysis and curve fitting to the Michaelis-Menten equation. All experiments and controls were performed in triplicate.
In silico analysis: All in silico experiments were performed using YASARA modelling suite (YASARA structure v.18).
AlphaFold Colab: For this work, AlphaFold Colab notebook was utilized for the prediction of protein's 3D structure. [24] The relaxation option was also enabled for the final prediction. The quality of the final structure was evaluated visually, checking for the proper orientation of the catalytic residues, as well as with the AlphaFold's platform results.
Molecular docking: Molecular docking was performed using the VINA algorithm [27] and the AMBER03 force field [28] at 50°C in vacuo. The simulation cell was 7 Å towards each axis direction from the Oγ atom of the catalytic serine residue. The protein's side chains were kept flexible during the docking procedure, creating a total of five protein ensembles with different side-chain rotamers. For each ensemble, 25 docking runs were executed, and the resulting ligand structures clustered when the ligand RMSD was lower than 5 Å. Prior to docking, energy minimization was performed for both the protein structure and the ligands using YASARA2 force field. All figures were prepared by using PyMOL v 0.99. The parameters that were used to identify the catalytically active docked conformations of each ligand were: 1) stability of the hydrogen bonding interactions with the oxyanion hole, 2) the distance between the catalytic serine Oγ atom and the substrate Cα to be less than 4 Å, 3) low free energies of binding. Molecular dynamics simulations: MD simulations were used to analyse and further confirm the binding poses that occurred from molecular docking. The boundaries of the simulation cell were set to periodic to avoid surface tension effects. The TIP3P solvation model was used for the simulation of explicit water molecules and the YASARA2 force field for the simulation of the complexes. [29] The setup of this work included the optimization of the hydrogen bonding network [30] to increase the protein's stability, as well as a pK a prediction so that the protonation states of the protein's residues would represent the model at the chosen pH of 7.2. The cell charge was neutralized by adding NaCl ions with a physiological concentration of 0.9 %. Afterwards, the system was subjected to energy minimization with the steepest descent and simulated annealing algorithms in order to remove clashes. The simulation was then run for a total duration of 50 ns and the time step was set to 1.25 fs for bonded interactions and 2.5 fs for nonbonded interactions at a temperature of 323.15 K and constant pressure of 1 atm (NPT ensemble) using the Particle Mesh Ewald algorithm for the interactions of long-range electrostatic forces. [31] Finally, the coordinates of all the atomic positions were stored every 20 ps and the recurring snapshots were plotted for the trajectory analysis.