Homology modeling and docking studies of a Δ9-fatty acid desaturase from a Cold-tolerant Pseudomonas sp. AMS8

Templates selection

BLASTP of the target protein was performed at the NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp) which showed 24 and 23% similarities to human integral membrane stearoyl-CoA desaturase (PDB ID: 4ZYO) and mouse stearoyl-coA desaturase (PDB ID: 4YMK), respectively. However, the mouse stearoyl-coA desaturase (PDB ID: 4YMK) was chosen as a template (based on its high resolution) to predict the three dimensional structure of the Δ9-fatty acid desaturase.

Structure prediction

The primary sequences of human (Uniprot ID: O00767) and mouse (Uniprot ID: P13516) desaturases were obtained from Uniprot protein databases. Moreover, the protein sequences of Δ9-fatty acid desaturases from several other Pseudomonas species were obtained from the GenBank. The transmembrane (TM) spanning region of Pseudomonas sp. AMS8 Δ9-fatty acid desaturase was predicted using a TM domain topology prediction program, CCTOP (http://cctop.enzim.ttk.mta.hu/). The CCTOP predicted the TM domains of the protein sequence based on the consensus of ten (10) different methods (Dobson, Reményi & Tusnády, 2015a; Dobson, Reményi & Tusnády, 2015b). The TM domains were modelled together with the remaining amino acid residues toward the C-terminus of the protein using MODELLER version 9.18 (Webb & Sali, 2014). Alignment input used in the MODELLER was derived from the pairwise alignments of both the template and model primary sequences using membrane proteins alignment tool (Stamm et al., 2013) whereas the secondary structure was predicted using PSIPRED tool (Buchan et al., 2013) and used as a guide to modelling the C-terminal domain.

Energy minimization and quality verification

The protein model generated by the MODELLER was ranked and scored using discrete optimised protein energy (DOPE) score. The top three models out of 50 models with the lowest DOPE scores were chosen and assessed using ERRAT and RAMPAGE server. The best model was selected for energy minimization to remove clashes between side chains using GROMACS and subsequently used in docking and molecular dynamics simulations. Further structural assessment was performed by simulating the homology models inside a membrane bilayer.

Active site prediction

Prior to docking simulations, the binding site for the Pseudomonas sp. AMS8 model was predicted using COACH (Buchan et al., 2013). The predicted active site was found in the vicinity of histidine rich region, which served as an already established potential binding site of the enzyme (this site was used for a targeted docking). Blind docking was also performed to bind palmitate on both the template and model structures.

Docking studies

Three dimensional (3D) structure of palmitic acid was obtained from the Pubchem (https://pubchem.ncbi.nlm.nih.gov). The energy minimized model of the Δ9-fatty acid desaturase and the palmitic acid (ligand) were prepared and used for molecular docking of the substrate onto the target proteins using Autodock tool Version 4.2 (Trott & Olson, 2010). Blind docking of the palmitate onto the modelled structure was performed using a pre-set simulation grid box size of 126 × 126 × 126 Å along the X, Y and Z axes and centred at 39.946, 40.191, 45.879 whereas the targeted docking grid box size was set to 70 × 70 × 60 Å dimension and centred at 43.946, 40.191, 33.879 of X, Y and Z coordinate, respectively. The docking simulations were performed for 100 runs using Lamarckian Genetic Algorithm (LGA). The results were evaluated using RMSD values, ligand-protein interactions, binding energy (ΔG_bind) as well as a number of conformations existed in a populated cluster. The charge values of metal ions have been indicated to play a crucial role in predicting correct docking simulations. Previous report showed that different charge values of magnesium ions were verified and found a charge of +2 was too high, causing it to diminish due to partial binding of the adjacent groups (Chen et al., 2007). A docking study conducted on dioxyhypusine hydroylase (DOHH) with its substrate in the presence of Fe-Fe confirmed that Fe charge of +1 produced substrate binding mode which sufficiently agreed with the experimental data (Han et al., 2015). Thus, in this study, iron charges of 0 (default charge in Autodock), +1 and +2 were tested in the docking simulations for comparison. The ligand-protein interaction was visualized using Pymol (Trott & Olson, 2010) and VMD (Humphrey, Dalke & Schulten, 1996).

MD simulations

The model structure of the Δ9-fatty acid desaturase was simulated in an embedded 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer. The protein-bilayer system was constructed using CHARMM-GUI Membrane builder (Jo et al., 2008). The atomistic MD simulations were performed using CHARMM36 force field (Huang & MacKerell, 2013) in GROMACS Version 5 (Pronk et al., 2013) within an integration time step of 20 fs. The simulation temperature was kept constant at 310 K by coupling the system to a heat bath using a Nose–Hoover thermostat with τ_T = 1 ps. Pressure was maintained at 1 atm using a Parinello-Rahman barostat and semiisotropic pressure with τ_P = 5 ps and a compressibility of 4.5e⁻⁵ bar ⁻¹. Long-range electrostatics was treated using particle mesh Ewald method with a cutoff of 12 Å. The 12 Å cutoff distance was used for van der Waals interactions. The systems were equilibrated for 1 ns restraining the C α atoms, followed by production runs of 50 ns each in triplicates. The data was analysed using GROMACS tools and VMD.

Sequence of the Δ9-fatty acid desaturase protein and template identification

The Δ9-fatty acid desaturase was isolated from a cold-tolerant Pseudomonas sp. AMS8 and functionally expressed in Escherichia coli as confirmed by GCMS analysis which showed an active enzyme capable of increasing the overall palmitoleic acid content of the recombinant E. coli. Based on the GCMS analysis, a profound increase of the amount of palmitoleic acid from 10.5 to 21% was observed at 20 °C (Garba et al., 2016a). The protein had a molecular weight of 45 kDa and 394 amino acids which was already deposited at NCBI (accession number: AMX81567). Multiple sequences alignments of the template, the cold-tolerant Pseudomonas sp. AMS8 Δ9-fatty acid desaturase and sequences from several other Δ9-fatty acid desaturases have revealed the three conserved-histidine boxes common to all membrane-bound desaturases in bacteria (Garba et al., 2016b; Li et al., 2009), fungi (Chen et al., 2013) and animals (Bai et al., 2015) (Fig. 1). The human integral membrane stearoyl-CoA desaturase (PDB ID: 4ZYO) and mouse stearoyl-coA desaturase (PDB ID: 4YMK) have been solved to a resolution of 3.25 Å and 2.6 Å, respectively. The two structures share moderate sequence identities of 24 and 23% with the Pseudomonas sp. AMS8, respectively. However, the mouse stearoyl-coA desaturase was chosen as the template based on its higher resolution.

Model of the Δ9-fatty acid desaturase

To correctly model the TM domain of a membrane protein, it is necessary to appropriately predict its TM spanning region. CCTOP used 11 TM prediction programs including some of the best TM domain predictor such as TMHMM and HMMTOP to predict the TM domain of the target protein. Most of the prediction programs predicted that the target protein has three TM spanning regions as detailed in Fig. 2A, which gave a consensus domains of TM1 (Leu13–Leu33), TM2 (Leu135–Ile159) and TM3 (Met162–Tyr181) of 20, 25 and 20 amino acid residues, respectively. However, the TM2 and TM3 were not aligned at the TM domains of the template. Thus, some manual adjustment of the TM2 and TM3 was performed to prepare the alignment input for the MODELLER. As the three dimensional (3D) structure of the template had four TM domains, only three TM3 domains were considered for modelling the Pseudomonas sp. AMS8 Δ9-fatty acid desaturase (predicted to have only three TM domains) using MODELLER (Fig. 2B). The structure with the lowest DOPE score was assessed and improved after energy minimization and subsequently used for further analyses.

Quality verification of the predicted structure

Quality of protein models are verified using various programmes such as ERRAT and Ramachandran plot which are freely available online servers (Colovos & Yeates, 1993). In this study, the predicted structure that has gone through energy minimization was verified for correctness using the ERRAT and Ramachandran plot and labelled as AMX8-em (shown in the supplementary files). The ERRAT programme showed an overall quality value of 65.021% for this structure (Fig. 3). In general, at lower resolution, more than 92% of the surveyed structures had more than 80% of their residues outside the 95% exclusion zone (Colovos & Yeates, 1993). RAMPAGE programme is used to check the overall sterio-chemical quality, local and residue-by-residue reliability usually shown on a Ramachandran plot. The programme shows the sterio-chemistry of the main-chain torsion angles Phi, Psi (φ, ψ) of a good protein model. The Ramachandran plot displays the polypeptide chain of a protein structure using the φ, ψ angles pair (Laskowski, Moss & Thornton, 1993; Mahgoub & Bolad, 2013; Ramachandran, Ramakrishnan & Sasisekharan, 1963). Figure 4 and Table 1 indicate that up to 91.6% of the residues fall within the most favoured regions, 6.8% in the allowed regions whereas only 1.6% residues are in the outlier regions, further confirming that the predicted model is of good quality.

Catalytic site of the predicted structure

Membrane-bound desaturases share an exceptional structural resemblance and a wide range of functionality. Three conserved-histidine boxes that are common to all classes of these enzymes function in binding two irons at the catalytic centre. The structural similarity has given an insight into their structure-function relationships (Meesapyodsuk et al., 2007). The predicted structure of Pseudomonas sp. AMS8 Δ9-fatty acid desaturase contains the three conserved-histidine boxes consisting of eight histidine residues at positions 1 (His34, His39), 2 (His71, His74, His75) and 3 (His206, His209, His210) from N to C-terminus of the enzyme shown in Fig. 5A and analysed in Fig. 5B. The conserved-histidine motifs are consistent with those observed during the multiple sequences alignments of the target sequence with the template corresponding to the already established catalytic centre of membrane-bound desaturases. The role of the eight histidine residues in the conserved histidine-rich motifs has been demonstrated through site-directed mutagenesis of rat stearoyl-CoA Δ⁹-desaturase whereas those residues flanking the conserved region have critical catalytic properties in plant FAD2 desaturases and related enzymes (Broadwater, Whittle & Shanklin, 2002; Meesapyodsuk et al., 2007; Shanklin, Whittle & Fox, 1994).

Docking studies

The membrane-bound Δ9-fatty acid desaturase uses activated oxygen molecule to create double bond between C-H bonds of saturated substrates. The enzyme particularly introduces a double bond at Δ9-position of saturated palmitic and stearic acids to produce palmitoleic and oleic acids, respectively, serving as the fundamental substrates for phospholipids construction and other complex lipid molecules (Castro et al., 2011).

To investigate substrate specificity of Pseudomonas sp. AMS8 Δ9-fatty acid desaturase, docking studies of palmitic acid onto the modeled structure and the template were performed using Autodock software. Blind docking of palmitate and the template was first performed which showed that the palmitate was docked on the template at a site different from the vicinity of the template catalytic site observed for its native ligand. Similarly, for the model structure of desaturase from Pseudomonas sp. AMS8, the docked conformation with lowest docking energy formed close contact with Thr4, Trp167, Val171, Gly170, Leu175, Ala63, Cys96, Gly174, Tyr95, Met141, Ile144, and Ile140 outside the catalytic site. This is expected as the COACH predicted multiple binding sites on the protein.

It is known that the potential catalytic and binding sites for palmitate are close to the His conserved motif. Therefore, specific docking was performed with grid box which covers the histidine residues of the motifs. The docking simulation which produced the lowest binding energy, −6.8 kcal/mol is depicted in Fig. 6A. A number of residues were found in close contact with the docked palmitic acid namely, Ile26, Tyr95, Val179, Gly180, Pro64, Glu203, His34, His206, His71, Arg182, Thr85, Lys98 and His177 (Fig. 6B). Interestingly, among the binding residues are His34 and His71 and His206 from the first, second, and third conserved histidine motif of the enzyme, respectively. The ligand formed two hydrogen bonds with Lys98 and His177. These suggest that the docked substrate was very close to the enzyme catalytic site and the conserved-histidine residues holding the metal ions of membrane-bound desaturases which are known to play a key role for the enzymes catalytic activity (Shanklin et al., 2009).

To ascertain the role of the metal irons on the membrane-bound fatty acid desaturases, docking simulations were performed on both the model (Fig. 7) and template (Fig. 8) structures in the presence and absence of metal irons at their catalytic sites. Using the Autodock, different docking scores were predicted and summarized in Table 2. Lowest docking energy (−6.81) was observed when the palmitate was docked onto the model in the presence of metal irons, which suggests a more favourable configuration (Fig. 9). However, the control docking simulations performed on the template with the palmitate in the presence and absence of the iron metals showed much higher docking energies of −6.21 and 5.82, respectively, which indicate that the palmitate was unfavourably bound to the mouse desaturase enzyme catalytic site even in the presence of the metal irons.

Simulation of the predicted model of the Δ9-fatty acid desaturase in membrane

In order to further assess the Δ9-fatty acid desaturase model, the protein was embedded within a POPC bilayer and simulated for 50 ns. The approximate location of the bilayer was predicted based on the position of bilayer of the template structure. During equilibration, the protein movement was restrained for 1 ns to allow lipid molecules to equilibrate around the protein. The initial structure of the protein inside the POPC bilayer and its final structure at the end of the MD simulation are shown in Fig. 10.

In the membrane, the protein root mean square deviation (RMSD) was calculated to check the overall structure stability (Fig. 11). The RMSD underwent major changes in the first 5 ns and became more stable after 20 ns. To ascertain the real cause of the significant changes observed during the MD simulation, a more fine-grained analysis of the RMSD of the entire protein and its different regions (specifically the TM domains 1 to 3, and the C-terminal region following the last residue of TM3 (182 to 252) of the model structure) along the simulation was carried out. The RMSD established an overall stable protein, which is well equilibrated in the POPC bilayer after the 20 ns. The pdb coordinates of the protein analysed at 25, 30 and 50 ns were extracted from the MD trajectory and performed ERRAT analyses which showed overall quality scores of 90.81, 93.33 and 84.74%, respectively (Fig. 12).