Desaturase specificity is controlled by the physicochemical properties of a single amino acid residue in the substrate binding tunnel.

Membrane fatty acyl desaturases (mFAD) are ubiquitous enzymes in eukaryotes. They introduce double bonds into fatty acids (FAs), producing structurally diverse unsaturated FAs which serve as membrane lipid components or precursors of signaling molecules. The mechanisms controlling enzymatic specificity and selectivity of desaturation are, however, poorly understood. We found that the physicochemical properties, particularly side chain volume, of a single amino acid (aa) residue in insect mFADs (Lepidoptera: Bombyx mori and Manduca sexta) control the desaturation products. Molecular dynamics simulations of systems comprising wild-type or mutant mFADs with fatty acyl-CoA substrates revealed that the single aa substitution likely directs the outcome of the desaturation reaction by modulating the distance between substrate fatty acyl carbon atoms and active center metal ions. These findings, as well as our methodology combining mFAD mutational screening with molecular dynamics simulations, will facilitate prediction of desaturation products and facilitate engineering of mFADs for biotechnological applications.

Introduction

Membrane fatty acyl desaturases (mFADs) are present in all eukaryotic organisms and some prokaryotes. They catalyze a highly energetically demanding oxygen-dependent extraction of hydrogen atoms from fatty acyl hydrocarbon chains [1]. Individual members of the mFAD gene family exhibit rather narrow substrate specificity with respect to the fatty acyl head group, chain length, and pre-existing double bonds. mFAD homologs also show distinct regio- and stereospecificities with respect to the position and configuration of the double bond(s) they introduce, respectively [2], [3]. The vast repertoire of mFAD specificities results in an immense number of unsaturated fatty acyls (UFAs) with diverse physical and chemical properties that underlie their biological roles as cell membrane lipid components or precursors of signaling molecules including prostaglandins, leukotrienes, and insect pheromones [2], [3], [4].

Despite the fundamental role of mFADs in determining the composition of UFA pools within organisms, the structural determinants of mFAD specificities remain poorly understood, chiefly due to the difficulties of experimentally determining membrane protein structures. Two recent breakthrough structures of mammalian mFADs indicate that the acyl chain of the substrate is enclosed in a sharply kinked hydrophobic binding tunnel, which presumably constrains the substrate chain length and position with respect to the mFAD active center. This substrate binding tunnel likely provides a basis for the specificity of the desaturation reaction. The active center, which activates molecular oxygen necessary for the desaturation reaction, is formed by a di-iron cluster coordinated by His residues [5], [6].

To uncover the structural and sequential motifs that govern the outcome of desaturation reactions, more than a hundred mFAD genes and mutants have been cloned and functionally characterized in a yeast expression system [5], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20]. These studies identified an extensive collection of mFAD mutations that influence mFAD specificity, yet insights into universal mechanisms of desaturation specificity remain limited. Ideally, a search for mFAD specificity determinants would 1) use an experimentally determined mFAD structure [5] or homology-based model [13], [15], [16], [17], 2) cross-validate the effect of the identified mutations in more distantly related mFADs, and 3) systematically screen the effect of a range of amino acid (aa) residues at the candidate specificity-determining aa sites. Such a systematic approach has not yet been applied to study mFAD structure–function relationships.

Homology-based modeling of mFADs using either of the two available mFAD structures [5], [6] represents a feasible and increasingly accurate alternative to experimental structure determination [21], [22]. Previously, mapping aa residues and regions critical for mFAD specificity onto mFAD structural models localized the specificity-determining regions to the vicinity of the substrate-binding channel opening [15], distal end of the substrate-binding channel [5], [13], kink of the substrate-binding channel [17], and vicinity of the iron atoms in the active center [15], [16].

To shed light on the mechanisms of specificity determination in mFADs, we used as a model system two pheromone-producing mFADs, MsexD2 and MsexD3, from the tobacco hornworm moth Manduca sexta. This pair of mFADs with high sequence similarity exhibits distinct desaturase selectivities and specificities with C14 and C16 fatty acyl-CoA substrate resulting in production of a spectrum of mono-, di- and tri-unsaturated fatty acyls (Fig. 1) [17]. Notably, reciprocal swapping of amino acid residue Ala224/Ile224 within the predicted substrate-binding tunnel leads to exchange of the desaturase specificities and selectivities of MsexD2 and MsexD3 [17]. However, the mechanism by which a single amino acid residue modulates the desaturase substrate preference and the position of introduced double bond(s) remains elusive. Here we combine 1) experimental characterization of mFAD mutants constructed by systematic replacement of a single aa position with a panel of residues, 2) assessment of the mutation effect in a more distantly related mFAD from the silkworm moth Bombyx mori, and 3) molecular dynamics (MD) simulations of homology-based mFAD-substrate complexes. Our results provide mechanistic insights into how steric and hydrophobic effects of aa residues in the mFAD substrate-binding channel influence the position of the substrate with respect to the active center and thus govern the outcome of the desaturation reaction.

Fig. 1

Overview of desaturation reactions catalyzed by M. sexta and B. mori wild-type mFADs. The products of desaturation of tetradecanoyl-CoA (A) and hexadecanoyl-CoA (B) substrates via wild-type mFADs from M. sexta and B. mori as determined in yeast expression system. The position of introduced double bonds and the enzymes catalyzing these steps are shown along the reaction arrows. Conjugation step refers to synthesis of conjugated double bond system from monounsaturated precursor via 1,4-dehydrogenation which is formally accompanied by substrate double bond rearrangement [51]. Based on references [17], [24] and this study.

Materials and methods

Desaturase cloning and site-directed mutagenesis

The empty expression plasmid was provided by Holz et al. [23]. The pYEXTHS-BN expression plasmids bearing MsexD2, MsexD3, MsexD2-Ile224, and MsexD3-Ala224 coding regions with N-terminal His6tags under the control of the Cu+2-inducible promoter CUP1 were previously constructed in our laboratory [17]. Coding regions of wt B. mori desaturase BmorD1 [24] and its single-aa mutant BmorD1-Ile227 were optimized for yeast codon usage, custom synthesized (GenScript), and subcloned into the pYEXTHS-BN vector using the primers BmorD1(opt)_BamHI_F and BmorD1(opt)_EcoRI_R (Table S1). The corresponding GenBank accession numbers are MG450705 and MG450706. Plasmids bearing single point mutations coding for Gly224, Thr224, Val224, and Phe224 were constructed by site-directed mutagenesis with Phusion HF DNA polymerase (New England BioLabs) and mutagenesis primers (for primer list, see Table S1) to introduce single point mutations into pYEXTHS-BN vectors bearing either the MsexD2 or MsexD3 ORF. The cycling conditions comprised initial denaturation at 98 °C (2 min); 30 amplification cycles consisting of the following steps: 98 °C (10 s), 58 °C (20 s) and 72 °C (6 min); and a final extension step at 72 °C for 10 min. The PCR reaction product was treated with DpnI endonuclease for 1 h at 37 °C, and 1 µl of the reaction mixture was transformed into competent E. coli. All vectors were verified by sequencing.

Functional expression in yeast

The pYEXTHS-BN constructs and an empty pYEXTHS-BN vector (control) were transformed into a desaturase- and elongase-deficient double mutant Saccharomyces cerevisiae strain (MATa elo1:HIS3 ole1:LEU2 ade2 his3 leu2 ura3; further referred to as elo1Δole1Δ) [25] using the S.c. EasyComp Transformation Kit (Invitrogen) according to the manufacturer’s instructions. Transformed colonies were selected on agar plates with synthetic complete media lacking uracil [SC-U: 0.67% yeast nitrogen base without amino acids, 2% glucose, 2% agarose, supplemented with Brent supplement mix without uracil (ForMedium), according to the manufacturer’s instructions] supplemented with a mixture of unsaturated fatty acids [palmitoleic and oleic acid, each 0.05 mM, solubilized by 0.25% tergitol (type NP-40, Sigma-Aldrich)] to maintain the growth of the Z9-UFA deficient elo1Δole1Δ yeast strain.

Yeast strains were cultivated in 20 ml liquid SC-U medium supplemented with palmitoleic and oleic acid (each 0.05 mM), 0.25% tergitol, and 0.5 mM CuSO4 to induce heterologous expression of mFADs under the control of CUP1 promoter. Where indicated, yeast strains were cultivated with methyl Z11-hexadecenoate (Z11-16:1) supplemented in the induction medium to a final concentration of 0.1 mM or methyl E10E12-hexadecadienoate (E10E12-16:2) supplemented to a final concentration of 0.25 mM. Fatty acids 14:0 and 16:0, which served as substrates for mFAD E/Z11-desaturase activities, were naturally abundant in yeast cells and were not supplemented into the cultivation media, with the exception of cultivations of BmorD1-transformed yeast cultures for analysis of C14 products, which were supplemented with 0.25 mM methyl tetradecanoate. The yeast liquid cultures were cultivated for 3 days until they reached the stationary phase and the yeast cell biomass was harvested by centrifugation as previously described [26]. Heterologous mFAD expression was monitored by western blotting with anti-His6tag antibody [27]. All mFADs migrated on SDS-PAGE as apparently lower molecular weight proteins compared to their calculated protein molecular weights of 37.2 kDa and 39.6 kDa, respectively, presumably due to anomalous migration of transmembrane proteins [28]. Total cellular lipids were transesterified to fatty acid methyl esters (FAMEs) and extracted as described by Matoušková et al. [29].

GC/MS analysis

Conditions for GC/MS analysis of extracts from yeasts expressing MsexD2 and MsexD3 variants, and from yeasts expressing BmorD1 supplemented with methyl tetradecanoate were as follows: carrier gas, He, 1 ml/min; splitless injection, 1 μl; injector temperature, 220 °C; thermal gradient, 40 °C for 2 min, 40 °C to 140 °C at 70 °C/min, 140 °C to 240 °C at 3 °C/min, 240 °C to 280 °C at 20 °C/min, and 280 °C for 5 min. The temperature program was terminated at 230 °C and held for 15 min when the DB-WAX column was used. The extracts from yeasts expressing BmorD1 (except for those from methyl tetradecanoate supplementation) were analyzed as follows: carrier gas, He, 1 ml/min; 1 μl split injection 5:1; injector temperature, 250 °C; DB–5 ms column; thermal gradient, 140 °C for 1 min, 140 °C to 230 °C at 3 °C/min, 230 °C to 320 °C at 20 °C/min, and 320 °C for 3 min.

The relative FAME quantity of a particular FAME was calculated as the total ion current (TIC) area of its peak divided by the sum of TIC areas of all major FAME peaks in the chromatogram; the resulting value was then multiplied by 100%. The differences between mean values of relative FAME production at a significance level of 0.05 were calculated using ANOVA with a “post-hoc” Tukey’s test using R language and agricolae v1.2–6 R-package [30].

Desaturase sequence analysis

Desaturase sequences retrieved from NCBI GenBank were aligned using the Muscle algorithm [31] and formatted using JalView [32].

Structure modeling, molecular dynamics simulations and quantum-chemical calculations

The initial structures of wt MsexD2 and MsexD3 were generated using the homology modeling module in YASARA [21], with default parameters and charges corresponding to pH 7 and using the structure of mammalian mFAD with PDB ID code 4YMK as a template [5], [17]. Wt MsexD2, corresponding mutants, and ions were described with the Amber99SB force field [33], while parameters for fatty acyl-CoA substrates were based on the Slipids force field [34], [35], [36]. All missing bonding and nonbonding parameters of fatty acyl-CoA substrates in the existing Slipids force field were updated with compatible CHARMM36 parameters, while atomic charges were calculated by the Merz-Singh-Kollman scheme [37] consisting of B3LYP/6-31G(d) geometry optimization and a subsequent single point ESP charge calculation using the B3LYP/cc-pVTZ method. All systems were solvated in a cubic box of dimensions 10.0 nm × 10.0 nm × 10.0 nm with approximately 31,500 TIP3p water molecules [38]. The total charge of the MsexD2 protein was +9, whereas the charge of the fatty acyl-CoA substrates was −4. Electroneutrality of the system was maintained by adding 5 chloride counterions. Two zinc atoms in the MsexD2 protein were constrained using umbrella potential with a force constant of 1000 kJ mol−1 nm−2 to keep the active site of the protein stable. Zn atoms were not covalently bound to histidine residues and were allowed to freely accommodate in the active site. The stability of the active site was checked by monitoring its hydration and water coordination around histidines, which was kept roughly constant throughout MD simulations. At this stage, simulations with MsexD3 were unstable for most mutants and the substrate was not optimally fitted in the active site, and therefore further MD simulations were performed with the focus on wt MsexD2 and five MsexD2 mutants with single aa substitutions at position 224 (Thr224, Val224, Ile224, Gly224 and Phe224) together with two different fatty acyl-CoA substrates, 16:0 and E10E12-16:2.

All MD simulations were run for 100 ns after initial equilibration of at least 10 ns with employed periodic boundary conditions in three dimensions. Long range electrostatic interactions beyond the non-bonded cutoff of 1.0 nm were described with the particle-mesh Ewald procedure [39] using Fourier spacing of 0.12 nm. The real space Coulomb interactions were cut off at 1.0 nm and van der Waals interactions were cut off at 1.4 nm. Pressure of 1 bar was maintained using the Parrinello–Rahman algorithm [40] with a coupling constant of 10 ps−1. The temperature of 310 K was controlled with the Nose–Hoover thermostat with a coupling constant of 0.5 ps−1 [41] independently for the protein with the fatty acyl substrate and water sub-systems. All bonds within protein and the fatty acyl-CoA substrate were constrained using the LINCS algorithm [42]. O–H bonds in water were kept constant using the SETTLE method [43]. Equations of motion were integrated using the leap-frog algorithm with a time step of 2 fs. All MD simulations were performed using the GROMACS program package, version 4.5.4 [44]. The number of molecular contacts between substrate carbon atoms and two Zn atoms in the active site below 0.7 nm, average distances, time evolution of distances, and standard deviations of average distances (substrate dynamics) were calculated for the duration of the MD simulations.

Results and discussion

aa224 in MsexD2 and MsexD3 fine tunes the desaturase activities and specificities

Our previous research highlighted aa224 residue as critical determinant of mFAD specificity in MsexD2 and MsexD3 by showing that reciprocal exchange of Ala224 and Ile224 residues which are present in the wild-type (wt) MsexD2 and MsexD3, respectively, leads to exchange of their desaturase specificities [17]. To uncover the mechanistic details of the specificity switch, here we first constructed a set of MsexD2 and MsexD3 mutants with glycine, alanine, valine, isoleucine, phenylalanine, or threonine at position 224 and expressed them in a yeast expression system. Western blotting of yeast lysates showed generally higher expression levels of MsexD3 and its mutants than that of MsexD2 and its mutants (Fig. 2). The expression levels among all wt and mutant MsexD2 and MsexD3, respectively, were comparable (Fig. 2), indicating that the aa224 mutations did not systematically or substantially influence mFAD protein levels in yeast cells. The relative production of UFAs by heterologously expressed mFADs was quantified by GC/MS analyses of fatty acyl methyl esters (FAMEs) transesterified from yeast total cellular lipids. The relative FAME quantities were then used to compare the approximate mFAD activity levels of the individual mutants and wt.

Fig. 2

Detection of heterologously expressed mFADs in yeast lysates by western blotting with anti-His6tag antibody. Yeast lysates from 50 μl of yeast culture were loaded in each lane. Anti-His6tag antibody was used to visualize the His-tagged mFADs. The bars on the left indicate positions of protein standard bands of 32 and 46 kDa. The lysate of a yeast strain transformed with an empty expression vector is marked “E.”

Mutations at aa224 in MsexD2 and MsexD3 led to a range of effects, including gain of novel or loss of original desaturase specificities, as well as increases or decreases in desaturase activity relative to the wt mFADs (Fig. 3). Desaturase activity was generally higher in MsexD3 and its mutants, probably as a result of their higher expression levels (Fig. 2). MsexD2- and MsexD3-aa224 mutants exhibited strikingly similar trends in desaturase specificity and activity with respect to the aa224 side chain volume (Fig. 3), leading us to interpret the changes in mFAD enzymatic specificities in terms of the steric effect of the aa224 side chain. However, other physicochemical properties of the aa224 side chain, such as its hydrophobicity, also are likely to contribute to the observed shifts in UFA products.

Fig. 3

Production of unsaturated fatty acids by wild-type and mutant MsexD2, MsexD3, and BmorD1 desaturases. Relative abundances of monounsaturated UFAs E11-14:1 (A), Z11-14:1 (B), and Z11-16:1 (C); di-unsaturated UFAs E10E12-16:2 (D) and E10Z12-16:2 (E); and tri-unsaturated UFAs E10E12E14-16:3 (F) and E10E12Z14-16:3 (G) in lipid extracts of yeast strains heterologously expressing mFADs. The amino acid residues are listed in order of their increasing mean volume, as adopted from Harpaz et al.[50]. “Empty” represents the yeast strain transformed with an empty expression plasmid. The relative fatty acyl abundances values are the average of three yeast cultivation experiments ± SD, and the colored points indicate individual data points. The mean values which are significantly different (alpha = 0.05) in ANOVA with “post-hoc” Tukey's test are marked with different letters.

For MsexD2, MsexD3, and their mutants, Z11-desaturation of the C16 substrate and conjugation activity (i.e., production of conjugated E10E/Z12-16:2 by a combination of desaturation and shift in preexisting double bond position) followed the same trend and were favored by less bulky Ala224 and Thr224 residues (Fig. 3D,E). In contrast, bulkier Val224 and Ile224 favored E/Z14-desaturase activity with the E10E12-16:2 substrate. The least bulky and the bulkiest aa224 residues, Gly and Phe, respectively, resulted in an overall decrease in or loss of desaturase activity (Fig. 3C).

Z/E11-desaturation of a short substrate (C14) in both MsexD2, MsexD3, and their mutants was favored by bulkier aa224 residues (Fig. 3A, B). In contrast to desaturation of the C16 substrate, Z/E11 desaturation of the C14 substrate did not dramatically decrease for the bulkiest Phe residue. These results suggest that steric constraints for desaturation of shorter fatty acyl substrates, such as C14, might be relaxed compared to fatty acyl substrates with longer chains.

The role of aa224 as a desaturase specificity switch is conserved in B. mori mFAD

To evaluate whether the role of aa224 in determining desaturase specificity also is conserved in other mFADs, we mutated BmorD1, a pheromone-producing mFAD from the silk moth Bombyx mori [24], which shares 62.9% and 60.2% sequence identity with MsexD2 and MsexD3, respectively (Fig. S1). BmorD1 desaturase specificity closely parallels that of MsexD2 by exhibiting both Z11-desaturase and conjugase activity [24]. Despite these similarities, the gene tree of mFAD gene family does not support one-to-one orthology between MsexD2 and BmorD1 [17]. Instead, it is plausible that the conjugase activities of MsexD2 and BmorD1 evolved independently and convergently.

Thr227 in BmorD1 is homologous to Ala224 in MsexD2 and Ile224 in MsexD3 (Fig. S1). Replacing this residue with Ile significantly decreased production of Z11-16:1 compared to wt BmorD1 and abolished production of E10E12-16:2 and E10Z12-16:2 (Fig. 3C-E). Remarkably, this mutation also conferred novel E/Z14-desaturase activity onto the BmorD1-Ile227 mutant (Fig. 3F, G). The BmorD1-Ile227 mutant thus acquired a desaturase specificity profile close to that of MsexD2-Ile224 and MsexD3-Ile224 (wt). These results indicate that mFAD specificity is to a large extent controlled by conserved mechanisms and can be tuned by changes as subtle as a single aa substitution in the kink of the FA-binding tunnel. The presence of Ile224 alone, however, is apparently not sufficient to install E/Z14-desaturase specificity across all lepidopteran mFADs, as several moth mFADs with Ile at the position homologous to aa224 do not exhibit E/Z14-desaturase specificity (Fig. S2) [45], [46], [47], [48], [49]. In these mFADs, structural changes in other regions likely further modulate the desaturase specificity.

Mutations of aa224 modulate the intermolecular distances between substrate and active center

To uncover the mechanistic details of the experimentally observed desaturase specificity switch, we performed MD simulations of complexes of mFADs with fatty acyl substrates. We analyzed the approach of a reactive fatty acyl substrate bond to the mFAD active site metal ions and overall substrate dynamics (quantified as standard deviation of average distances between substrate carbon atoms and Zn atoms in the active site).

First, we generated homology models of MsexD2 and MsexD3 with substrates 16:0 and E10E12-16:2 inserted in the substrate binding tunnel, using mammalian mFAD structure as a template [5], [17]. Our models included active center Zn ions rather than the in vivo catalytically active iron ions [1] to comply with the experimentally determined mammalian mFAD structures [5], [6]. In the MsexD2 models, the C11-C12 atoms of 16:0 and the C14-C15 atoms of E10E12-16:2, where desaturation occurs, are positioned in the vicinity of two active center Zn ions. The fatty acyl chain is deeply inserted into the substrate binding tunnel and wrapped around the kink in the tunnel formed by aa224 (Fig. 4A, E). We performed MD simulations with both MsexD2 and MsexD3. However, the MD simulations of MsexD3 were unstable and we therefore proceeded with simulations of MsexD2 and its mutants. Since the panel of amino acid residues introduced at position aa224 had comparable effect on desaturase specificity in both MsexD2 and MsexD3 backgrounds, we assume that either of MsexD2 and MsexD3 are suitable for simulations of the desaturase specificity switch. Nevertheless, the instability of MsexD3 in MD simulations is surprising, given the high overall sequence identity of MsexD3 and MsexD2 [91% in the homologous 321-aa region (Fig. S1)]. The instability of MsexD3 model in MD simulations could be explained by increased constrain of the substrate imposed by bulkier Thr223 in MsexD3 (as compared to Ala223 in MsexD2). Assuming certain level of force field inaccuracy in the MD simulations, Thr223 could excessively constrain the substrate in the MD models and render the complex unstable. We have previously demonstrated that replacement of Ala223 by Thr223 in MsexD2 decreases its overall desaturase activity [17] which supports the possible role of Thr223 as a constraint in substrate binding.

Fig. 4

Molecular dynamics simulations of MsexD2 with 16:0-CoA (A, B, C) and E10E12-16:2-CoA (D, E, F) substrates in the substrate-binding tunnel. Snapshots from MD simulations of MsexD2 with 16:0-CoA (A) and E10E12-16:2-CoA (D) substrates show Ala224 in van der Waals representation; Zn atoms in the active site are shown in grey surrounded by His residues and water molecules. Atoms belonging to the single bond of the substrate where desaturation takes place are labeled, along with the corresponding distance to Zn atoms in the active site. The total number of contacts (less than 0.7 nm) between the Zn atoms in the enzyme active center and the C11-C12 atoms of 16:0-CoA (B) and the C14-C15 atoms of E10E12-16:2-CoA (E) are indicated. The average distances ± SD between Zn atoms and desaturated carbon atoms are plotted together with relative unsaturated fatty acyl production (C,F) reproduced from Fig. 3, with SD bars omitted for clarity. The average distances which are significantly different (alpha = 0.05) in ANOVA with “post-hoc” Tukey's test are marked with different letters.

MD simulations of MsexD2 with 16:0 substrate

Calculations indicated that MsexD2-Thr224 and wt MsexD2 (Ala at position 224) with the 16:0 substrate exhibits the shortest average distances between the active site Zn atoms of mFAD and the carbon atoms of the C11-C12 single bond, as well as the largest total number of contacts (Table S2, Fig. 4B,C). These MsexD2 variants also exhibited the highest Z11-desaturase activity toward the 16:0 substrate among the experimentally characterized mFADs (Fig. 3C). The next shortest average Zn/C11-C12 distances and numbers of Zn/C11-C12 contacts were calculated for MsexD2–Ile224 and MsexD2-Val224 (Fig. 4B,C). Notably, MD simulation statistics of MsexD2–Ile224 with 16:0 substrate are suggesting higher desaturase activity than detected in experiments. This disagreement could be caused by our inability to fully describe the desaturase specificity mechanism in terms of simple MD simulation statistics. In MsexD2-Gly224, the average distance between the active site and substrate carbon atoms was large (greater than 1 nm), in agreement with the lower Z11-desaturase activity of this mutant. The large average Zn/C11-C12 distances calculated for the MsexD2-Phe224 complex with 16:0 suggest that the bulky, hydrophobic Phe side chain hinders the approach of C11-C12 atoms to the active center (Table S2, Fig. 4C), explaining the loss of desaturation activity (Fig. 3C).

In summary, the numbers of contacts and average distances between the C11-C12 atoms of the 16:0 substrate and the Zn active center ions are in generally good agreement with the experimentally determined Z11-desaturase activity of wt MsexD2 and its mutants. We observed disagreement of MD simulation and experiments for MsexD2–Ile224 with 16:0 substrate which however does not invalidate our ability to identify optimal aa224 residues for Z11-desaturase activity using MD simulation statistics.

MD simulations of MsexD2 with 16:2 substrate

In MD simulations of wt and mutant MsexD2 with the E10E12-16:2 substrate, MsexD2-Val224 exhibited the shortest distance and the highest number of intermolecular contacts between the C14-C15 atoms of E10E12-16:2 and the Zn ions (Fig. 4E, F). These MD parameters are in agreement with the experimental finding that MsexD2-Val224 has the highest E/Z14-desaturase activity of the mFADs characterized. The next shortest distance between C14-C15 and Zn ions was observed for MsexD2-Ile224, which also gained E/Z14-desaturase specificity (Fig. 3F,G). MsexD2-Ala224 and MsexD2-Phe224 exhibited very large distances and no contacts between C14-C15 and Zn, in good agreement with the absence and very low E/Z14-desaturase activity observed for wt MsexD2 and MsexD2-Phe224, respectively (Table S2, Fig. 3F,G and Fig. 4E F). The relatively close approach of the C14-C15 substrate bond to the Zn ions and higher substrate dynamics (Table S3, Fig. S4C) may explain retention of low E/Z14-desaturase activity in MsexD2-Thr224 (Fig. 3F,G). MsexD2-Gly224 with the E10E12-16:2 substrate exhibited an intermediate average distance and similar number of C14-C15-Zn contacts as MsexD2-Ile224 (Fig. 3E,F and Fig. S4A, C, Table S3) but did not exhibit any detectable E/Z14-desaturase activity (Fig. 3F,G). Quantitative disagreement of MD simulations and experiments for MsexD2-Gly224 with E10E12-16:2 substrate however does not interfere with our ability to predict Val224 and Ile224, respectively, as optimal residues for E/Z14-desaturation of E10E12-16:2 substrate.

Together, the MD simulations revealed that the average distances between the active center metal ions and the substrate carbon atoms at which desaturation is predicted to occur, along with the number of contacts and substrate dynamics, are remarkably good predictors of mFAD specificities. Our experiments and MD simulations can be interpreted in terms of steric effects, as we identified distinct optimal aa224 side chain volumes for Z11-desaturase activity with the 16:0 substrate and E/Z14-desaturase activity with E10E12-16:2. For desaturation of the 16:0 substrate, the optimal aa224 side chain volume was 26–56 Å3 (Ala-Thr), while for the E10E12-16:2 substrate, the optimum was shifted towards more bulky Val and Ile residues with side chain volume ∼ 75–101 Å3 (see Harpaz et al. for amino acid side chain volume calculations [50]). The approach of substrate carbon atoms to the mFAD active center ions is likely also modulated by hydrophobic effects of the aa224 side chain. The aa224 residues that bring substrate most closely to the active center metal ions through steric and hydrophobic interactions exhibit the highest desaturase activity. Future studies using additional mFADs and expanded panels of amino acid residues will likely further refine the relative contributions of hydrophobic and steric effects to the approach of substrate carbon atoms to the active center.

Conclusion

In summary, we show that a combination of 1) candidate mutation screening with a panel of amino acid residues, 2) cross-validation of the mutation effect in distantly sequentially related mFADs, and 3) molecular dynamics simulations has the potential to uncover mechanistic details of desaturation specificity. Our results indicate that the side chain volume of a single amino acid residue in the mFAD binding tunnel controls the approach of substrate carbon atoms to the active center metal ions. This single residue thus directs the outcome of the mFAD-catalyzed desaturation. Our approach could help compensate for the scarcity of experimentally determined mFAD protein structures by predicting desaturase reaction outcomes from the mFAD protein sequence and ultimately enabling engineering of mFADs with desired enzymatic properties for biotechnological applications.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.