Molecular BioSystems View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
PAPER
Cite this: Mol. BioSyst., 2014, 10, 2942
View Journal | View Issue
Structural and functional investigation of zebrafish (Danio rerio) NOD1 leucine rich repeat domain and its interaction with iE-DAP† Jitendra Maharana,*a Bikash Ranjan Sahoo,b Aritra Bej,b Mahesh Chandra Patra,bc Budheswar Dehury,d Gopal Krushna Bhoi,b Santosh Kumar Lenka,b Jyoti Ranjan Sahoo,b Ajaya Kumar Routb and Bijay Kumar Beheraa Nucleotide binding and oligomerization domain 1 (NOD1), a cytoplasmic pattern recognition receptor (PRR) and is a key component for modulating innate immunity and signaling. It is highly specific to g-D-Glu-mDAP (iE-DAP), a cell wall component of Gram-negative and few Gram-positive bacteria. In the absence of the experimental structure of NOD1 leucine rich repeat (NOD1-LRR) domain, the NOD signaling cascade mediated through NOD1 and iE-DAP interaction is poorly understood. Herein, we modeled 3D structure of zebrafish NOD1-LRR (zNOD1-LRR) through a protein-threading approach and structural integrity of the model was assessed using molecular dynamics simulations. Molecular interaction analysis of iE-DAP and zNOD1-LRR, their complex stability and binding free energy studies were conducted to anticipate the ligand binding residues in zNOD1. Our study revealed that His775, Lys777, Asp803, Gly805, Trp807, Asn831, Ser833, Ile859 and Trp861 situated in the b-sheet region of zNOD1LRR could be involved in iE-DAP recognition, which correlates the earlier findings in human. Comparison
Received 3rd April 2014, Accepted 29th July 2014
of binding free energies of native and mutant zNOD1–iE-DAP complexes delineated His775, Lys777,
DOI: 10.1039/c4mb00212a
This study provides the first comprehensive description of biophysical and biochemical parameters
Asp803, Ser833 and Ile859 as the pivotal residues for energetic stability of NOD1 and iE-DAP interaction. responsible for NOD1 and iE-DAP interaction in zebrafish, which is expected to shed more light on
www.rsc.org/molecularbiosystems
NOD1 signaling and therapeutic applications in other organisms.
Introduction Innate immune system is the most ancestral and ubiquitous host defense mechanism to combat a wide range of pathogens. It is activated upon specific recognition of pathogen-associated molecular patterns (PAMPs) by pattern recognition receptors (PRRs).1,2 PRRs are widespread in extracellular, membrane and cytoplasmic regions, and are classified according to their PAMP specificity, function, localization and evolutionary origin. The nucleotide binding and oligomerization domain (NOD)-like
a
Biotechnology Laboratory, Central Inland Fisheries Research Institute, Barrackpore, Kolkata-700120, West Bengal, India. E-mail: [email protected]; Tel: +91 3325921190/91 b BIF-Centre, Department of Bioinformatics, Orissa University of Agriculture and Technology, Bhubaneswar-751003, Odisha, India c Animal Genomics Laboratory, Animal Biotechnology Centre, National Dairy Research Institute, Karnal-132001, Haryana, India d Department of Life Science and Bioinformatics, Assam University, Silchar-788011, Assam, India † Electronic supplementary information (ESI) available. See DOI: 10.1039/ c4mb00212a
2942 | Mol. BioSyst., 2014, 10, 2942--2953
receptors (NLRs) are some of the various reported PRRs, which belongs to the family of intracellular PRRs having canonical tripartite domain architecture. This includes an N-terminal effector binding domain (EBD), a central NACHT domain mediating self-regulation and oligomerization3 and a C-terminal leucine rich repeat (LRR) domain that recognizes PAMPs.4,5 The EBD domain is composed of either one or two caspase activating and recruitment domains (CARD),6 or a baculovirus inhibitor of apoptosis protein repeat (BIR),7 or a pyrin domain (PYD).8,9 Based on the diversity of the N-terminal EBD, NLRs are divided into five major subgroups, viz., NLRA (NLR family, acidic domain containing proteins), NLRB (NACHT, LRR and BIR-containing proteins), NLRC (NACHT, LRR and CARD-containing proteins), NLRP (NACHT, LRR and PYD-containing proteins) and NLRX (NLR family with no strong homology to the N-terminal domain of any other NLR subfamily member).10 The NLRC or NOD sub-group is constituted of five receptor proteins, NLRC1 (NOD1), NLRC2 (NOD2), NLRC3 (NOD3), NLRC4 (IPAF) and NLRC5 (NOD4), which sense various microbial derivatives and activate a complex downstream signaling cascade. Among these, NOD1 recognizes peptidoglycan (PGN)
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
derivatives, such as g-D-Glu-mDAP (iE-DAP), produced by all Gram-negative and specific Gram-positive bacteria.11–13 The NOD–PAMP interaction involves ATP recruitment which further leads to oligomerization of the NACHT domain, followed by CARD–CARD interaction between NOD and receptor-interacting serine/threonine-protein kinase 2 (RIP2), that engenders NF-kB activation and signal transduction.14 Previous studies showed that, NOD1 plays a significant role in innate immunity of higher eukaryotes.12,13,15 NOD1 has been recently characterized in lower eukaryotes like fish, frogs, snakes, etc. Notably, the role of the Nod1 gene in host protection has been studied in several bony fishes including zebrafish (Danio rerio),16 catfish (Ictalurus punctatus),17 grass carp (Ctenopharyngodon idella)18 and rohu (Labeo rohita).19,20 Over the past decade, with the advent of next-generation sequencing technologies, considerable advancements have been made in unraveling the mechanism of host–pathogen interactions. But, lack of structural information on NLR family members possesses a constraint to understand pathogen recognition and immune signaling pathways at the molecular level in higher-eukaryotes. In this scenario, zebrafish is recognized as an attractive model organism to study the in vivo importance of NLR proteins. In a recent study, Meeker et al.21 have shown that zebrafish shares notable resemblance with the human immune system. Furthermore, a study of Oehlers et al.16 has revealed that NLR protein family members, particularly NOD1 and NOD2 in zebrafish, are responsible for inflammatory bowel disease (IBD) predisposition. In this study, we have constructed the three dimensional (3D) model of zebrafish NOD1-LRR (zNOD1-LRR), using various molecular modeling techniques that includes multiple-template homology modeling, protein threading and ab initio modeling. We have established the low energy binding modes of iE-DAP and elucidated intermolecular interactions between iE-DAP and zNOD1-LRR using molecular docking, followed by molecular dynamics (MD) simulations. The binding free energy between zNOD1-LRR and iE-DAP was calculated using the Molecular Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) method. Taken together, this study for the first time demonstrates structural and dynamic properties of zNOD1-LRR in association with iE-DAP.
Materials and methods Domain analysis The protein sequence of zNOD1 was retrieved from the NCBI protein database (XP_002665106). The Conserved Domain Database (CDD) was used for analyzing the domain architecture of zNOD1.22 In addition to this, the zNOD1 protein sequence was also scanned in InterProScan,23 SMART24 and Pfam25 database to locate N-terminal CARD, central NACHT and C-terminal LRR domains. The LRR regions were manually identified by aligning the zNOD1-LRR domain with human, mouse and rohu. The multiple sequence alignment was executed in MAFFT26 and LRR regions were manually identified by locating ‘‘LxxLxLxxNxL’’
This journal is © The Royal Society of Chemistry 2014
Molecular BioSystems
motif (where L = leucine/isoleucine/valine/phenylalanine; x = any amino acid; N = asparagine/threonine/serine/cysteine). Template identification and molecular modeling To find suitable templates for the modeling of zNOD1-LRR, DELTA-BLAST27 search was performed against the protein data bank (PDB) (http://www.pdb.org/). Because of low sequence identity between target protein and PDB templates, the protein sequence was submitted to various protein threading servers, such as Genesilico Metaserver2,28 ModLink+29 and SPARKS-X,30 to identify suitable templates through fold recognition methods. Three best scoring templates resulting from fold recognition servers were considered for model construction using MODELLER 9.12.31 We have adopted both single-template and multipletemplate approaches for building zNOD1-LRR structures. The 3D models were also constructed using the automated modeling programs, SWISS-MODEL,32 BioSerf, (PS)2,33 M4T,34 EasyPred3D,35 Raptor-X,36 LOMETS37 and I-TASSER38 for comparison. The MODELLER-derived 3D models of zNOD1-LRR were evaluated based on discrete optimized protein energy (DOPE) scores. Three models with lowest DOPE scores were chosen for further refinement using GROMACS 4.5.5 simulation package.39 Molecular dynamics simulations MD simulations of the zNOD1-LRR models were carried out using Gromos96-53a6 force field and SPC/E water models in separate cubic boxes with a minimum distance of 15 Å between the protein surfaces and box edges. To neutralize the simulation systems, a physiological ionic strength (0.15 M NaCl) of counter ions was added. The atomic composition of the simulation systems is provided in Table S1 (ESI†). The final MD simulation of 20 ns was carried out for each system and the simulation parameters were followed from the previous report.40 The particle mesh Ewald (PME) summation method41 was used to treat the long-range electrostatic interactions. The non-bonded interactions were truncated after 12 Å of the cut-off radius. The temperature was coupled to 300 K using the Berendsen thermostat algorithm. A time step of 2 fs was used and the non-bonded neighbor list was updated heuristically. Conformations were saved at 2 ps time intervals for analysis. The average structure was extracted from the most stable portion of dynamics trajectory followed by energy minimization. The best model was chosen based on the stability of backbone atoms, individual residues, radius of gyration (Rg), secondary structure conservation and residual fluctuation. Model validation The stereochemical quality of the final energy-optimized model was verified using SAVeS (http://nihserver.mbi.ucla.edu/SAVES/), ProSA42 and ProQ.43 Bond length and bond angle analyses of the model was performed in MolProbity.44 Furthermore, the Z-score of hydrogen bond (H-bond) energy, packing defects, bump score, radius of gyration and deviation of Y angles of the refined model was tested in VADAR,45 GeNMR46 and PROSESS47 web servers. Time-dependent secondary structure analyses of the models were executed using visual molecular dynamics (VMD 1.9.1).48
Mol. BioSyst., 2014, 10, 2942--2953 | 2943
View Article Online
Molecular BioSystems
Paper
Two dimensional (2D) plots were generated using Grace 5.1.23 programs (http://plasma-gate.weizmann.ac.il/Grace/).
EMM = Eint + Ecoul + EvdW
Molecular docking of iE-DAP
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
The molecular mechanics interaction energy (EMM), is defined as (eqn (3)):
The 2D structure of iE-DAP (CID_45480617) was obtained from the PubChem database (http://pubchem.ncbi.nlm.nih.gov/). The Automated Topology Builder (ATB) server49 was used to generate the 3D coordinates of iE-DAP employing energy minimization with the Gromos96-53a6 force field. Molecular docking was carried out using AutoDock 4.250 to identify key ligand binding residues. The zNOD1-LRR model was assigned with Kollman charges and Geister partial charges were applied for iE-DAP using Autodock Tools 1.5. The binding site information of zNOD1-LRR was obtained from previously published literature.11,20 Initially, blind docking was performed to identify probable iE-DAP interaction sites. The docking was performed with different grid dimension around various LRR regions with a grid spacing of 0.375 Å. The initial population size was 300 and maximum number energy of evaluation was set to 2.5 107. The docked conformations were clustered with a 2 Å cut-off root mean square deviation (RMSD), and the resulting docked poses were ranked according to their binding energy scores. The best complexes from each docking simulation were subjected to MD simulations to optimize the ligand–receptor interaction.
where Eint represents bond, angle and torsion angle energies; Ecoul stands for electrostatic energy; and EvdW denotes van der Waals energy. The solvation free energy term (Gsol), is calculated using eqn (4): Gsol = Gpolar + Gnonpolar
A total of six docked zNOD1-LRR–iE-DAP complexes were selected on the basis of binding energy scores and the number of intermolecular H-bonds. Three independent MD simulations were performed on the selected docking complexes for better sampling of results. The ligand topology was prepared using ATB server.49 The procedure and parameters of MD simulations were the same as described in our previous work.40 Position restraints were applied to the backbone atoms of the modeled protein during equilibration simulation lasting 100 ps with both NVT and subsequent NPT conditions. The temperature was weakly coupled at 300 K via Berendsen method. The pressure was maintained at 1 bar through Parrinello–Rahman barostat algorithm in the NPT ensemble. After temperature and pressure equilibration, a total of 10 ns production MD was carried out for each complex. To amplify the performance and scaling of production recording, each simulation was truncated with three equal time intervals. For all six complexes, a total of 18 truncated MD (3 for each complex) was performed and the trajectories were analyzed. Binding free energy calculation A total of 500 snapshots were obtained from the MD trajectories to calculate the binding free energy51 by employing the MM/PBSA method52 as described in eqn (1): DGbinding = Gcomplex (GzNOD1 + GiE-DAP)
(1)
The individual energy terms for protein (zNOD1), ligand (iE-DAP) and complex are computed based on eqn (2): hGi = hEMMi + hGsoli ThSMMi
2944 | Mol. BioSyst., 2014, 10, 2942--2953
(2)
(4)
where Gpolar and Gnonpolar terms were calculated using APBS program.53 The polar term (Gpolar) was calculated by solving the non-linearized PB equation. The parameters for APBS calculation are as follows: grid spacing, 0.5 Å; temperature, 296 K; and salt concentration, 0.15 M. The nonpolar contribution, Gnonpolar, was computed using eqn (5): Gnonpolar = gSASA + b
(5)
where g = 0.0227 kJ mol1 Å2 and b = 0 kJ mol1.54 The dielectric boundary was defined using a probe radius of 1.4 Å. The standard errors were calculated using the following eqn (6): Standard error (SE) = s/ON
MD simulation of complexes
(3)
(6)
where s is the standard deviation and N is the number of structures used in the calculation.
Results and discussion Domain identification The domain analysis revealed that, zNOD1 sequence (940 amino acids (aa)) is comprised of a short N-terminal CARD (13–103 aa), a central NACHT domain (190–525 aa) and a C-terminal LRR domain (689–933 aa). Multiple sequence alignment of NOD1 sequences of zebrafish, human, mouse and rohu delineating the functional domains is shown in Fig. S1 (ESI†). The LRR domain in zNOD1 consisted of a series of covalently linked 21–29 aa residue motifs with the conserved pattern ‘‘LxxLxLxxNxL’’. Structurally, each LRR motif is comprised of one b-sheet with the pattern ‘xxLxLxx’ and one a-helix connected by a loop.55 Automatic prediction of LRR regions using SMART, Pfam and InterProScan advocated six LRR motifs in zNOD1. However, manual identification of LRR motifs from the aligned sequences of NOD1 from different species showed three additional LRR motifs (Fig. 1).20,40,56,57 Molecular modeling of zNOD1-LRR The search for templates for zNOD1-LRR in PDB database showed a very low percentage of sequence identities (i.e., o30%; twilight zone).58 So, the modeling exercises were performed employing fold recognition or threading methods.30 The resultant templates suggested by different threading programs are shown in Table S2 (ESI†). The 3D coordinates of zNOD1-LRR was assigned using crystal structures of porcine (PDB ID: 2BNH59),
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
Molecular BioSystems
Fig. 1 Multiple sequence alignment of human, mouse, rohu and zebrafish NOD1-LRR sequences using MAFFT. The leucine rich regions (LRRs) of zebrafish NOD1 are highlighted in boxes and the nine LRRs motifs are shown at the top of the alignment. The symbols ‘*’, ‘:’ and ‘.’ represents identical, conserved and semi-conserved substitutions, respectively.
mouse (3TSR) and human (1Z7X60) ribonuclease inhibitors. Models obtained from the automated molecular modeling servers were compared with those of MODELLER. To investigate the accuracy of built models, PSIPRED61 generated secondary structures of zNOD1-LRR were compared with the above mentioned template structures (Fig. S2, ESI†). After careful inspection, three top scored models were selected [model-1 (multiple templates); model-2 (Raptor-X); and model-3 (I-TASSER)] based on their DOPE scores (Fig. S3, ESI†). Analysis of MD trajectories The 20 ns MD trajectories of the selected zNOD1-LRR models were studied to gain insight into the structural characteristics. First, RMSD of backbone atoms was calculated as a function of simulation time. Both model-1 and 2 showed an average backbone deviation of 3.0 Å from its initial conformation; however, model-3 experienced an unusual backbone deviation (B6.5 Å) between 10 and 20 ns (Fig. 2A). The overall stability of zNOD1LRR was further confirmed by calculating Rg as a function of simulation time. The constant gyradius of 19.5 Å suggested a consistent shape and size of model-1 during simulation (Fig. 2B). In addition, the root mean square fluctuations (RMSF) of Ca atoms were calculated for each residue. The Ca atoms of
This journal is © The Royal Society of Chemistry 2014
Fig. 2 Evaluation of stability parameters of the modeled zNOD1-LRR in the three independent MD simulations (20 ns). (A) Root mean square deviations (RMSD), (B) radius of gyration (Rg) and (C) root mean square fluctuations (RMSFs). The color legends for different models are shown above the plot.
Mol. BioSyst., 2014, 10, 2942--2953 | 2945
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Molecular BioSystems
terminal amino acids fluctuated up to B10 Å in all models. The ligand binding domain C-terminal as studied earlier showed a prominent motion in model-1. LRR 6–8 regions in model-1 presented a distinct flipped conformation and indicated its engagement in ligand binding. On the other hand in model-2 and 3 these regions were comparatively stable (Fig. 2C). The C-terminal loop and helical regions of LRR domain were flexible during 20 ns simulation and suggested a structural rearrangement of functional domains in NOD1 upon PAMPs recognition.20 Secondary structure evolution from the trajectories showed that in model-1 all nine a-helices and nine b-strands retained their secondary structural properties throughout the MD simulation in aqueous solution (Fig. 3). Moreover, conformational dynamics analysis revealed no significant conformational alterations in model-1 during MD simulations as compared to model-2 and 3 (Fig. S4, ESI†). All together the MD simulation of top three models advocated model-1 as the most stable conformation and hence considered for further structural analysis. Structural analysis of the zNOD1-LRR model The optimized model obtained after MD simulation was investigated to delineate its key LRR domain and spatial conformation. The modeled zNOD1-LRR was displayed as a semicircle with a concave surface containing nine parallel b-strands, a convex face constituting nine a helices and short loops connecting these a/b-subunits. The LRRs within zNOD1-LRR are well conserved across different species as revealed from the sequence
Paper
Fig. 4 Schematic representation of the zNOD1-LRR model showing different LRR regions. The a-helices are shown in red, b-sheets in yellow and turns in green. Alignment of nine LRR motifs with the conserved pattern ‘‘LxxLxLxxNxL’’ is shown. The a-helices and b-sheets are displayed in red and blue charterers, respectively. The iE-DAP binding sites evidenced from previous reports are shown at the right bottom (transparent surface representations).
alignment (Fig. 4). Our study demonstrated that, the ligand binding active site of zNOD1 lies within LRRs 4–7 of the b-sheet region and corresponded to the previous reports in human (LRRs 5–7) by Girardin et al.11
Fig. 3 Secondary structures of three different models (A) model-1, (B) model-2 and (C) model-3 as a function of simulation time. Color segments given at the bottom of the figure represents the different secondary structural properties.
2946 | Mol. BioSyst., 2014, 10, 2942--2953
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Molecular BioSystems
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Model validation The stereo-chemical properties of the proposed zNOD1-LRR model were validated using various protein structure validation servers. The accuracy of dihedral angles (F/C) was gauzed by drawing the Ramachandran plot62 in PROCHECK integrated at SAVeS server. The zNOD1-LRR model showed that 77.90 and 98.60% of residues fell in the most favored and allowed regions, respectively, with an overall G-factor score of 0.04 (cut-off score 40.50) (Table 1). However, 1.40% of residues that includes Lys743, Gln864 and Ser868 fell within the disallowed regions. A close inspection of the modeled structures and previous reports revealed that the disallowed residues neither Table 1
Model validation report of zNOD1-LRR
Servers
zNOD1-LRR
PROCHECK
Verify3D ERRAT ProSA ProQ MolProbity VADAR
GeNMR PROSESS
Table 2
Most favored regions (%) Additionally allowed regions (%) Generously allowed regions (%) Disallowed regions (%) Overall G-factor Averaged 3D–1D score 40.2 Overall quality Z-Score LGscore MaxSub Cb deviations 40.25 Å Residues with bad bonds (%) Residues with bad angles (%) Standard deviation of w1 pooled Mean H-bond energy Generously allowed O angles (%) Packing defects (%) Percentage of 95% buried residues Ramachandran outside of most favored Bump score Radius gyration score Deviation of Y angles w1 score
77.90 18.90 1.80 1.40 0.04 99.59 75.10 7.09 6.33 0.54 0.00 0.00 0.00 2.77 0.70 0.26 0.21 0.54 2.56 0.34 1.20 1.23 1.20
involved in ligand binding nor constituted the structurally rigid regions. This indicated that the stereo-chemical properties of the model were within the acceptable range.63 The Verify3D score of the proposed model was 99.59%, which is above the cut off (Z80%) value and indicated that the quality of the proposed model is reasonably good.64 The ERRAT score provides accuracy of the non-bonded atomic contacts and our predicted model has a score of 75.10%, which is greater than the acceptable value (50%).65 The Z-score of our proposed model calculated using ProSA was in agreement with that of the experimental PDB structures of similar sizes. ProQ analysis indicated that the quality of zNOD1-LRR model was ‘extremely good’ based on their MaxSub and LGscores. Bond length and bond angle analysis of the proposed model using MolProbity revealed that none of the residues had bad side chain or main chain conflicts. The H-bond energy, packing defect, bump score, Rg and deviation of Y angles of the proposed model as analyzed by VADAR, GeNMR and PROSESS servers were within the cutoff range. A detailed validation report is summarized in Table 1. Analysis of docking results To identify the possible iE-DAP binding site in zNOD1-LRR, several molecular docking simulations were performed. Earlier studies have revealed that in NOD1 and NOD2, the b-sheet regions of the C-terminus (hot spot) appears to be critical for muropeptide recognition.11,40,56,66 In addition, a recent computational study on rohu NOD1 predicted C-terminal LRRs 8–9 as the iE-DAP binding site.20 Keeping all the evidences in support of the importance of C-terminal conserved LRR domain for iE-DAP recognition encouraged us to further investigate the binding modes between zNOD1-LRR and iE-DAP. In this study, we performed a total of 21 independent docking calculations to investigate the interaction of iE-DAP with zNOD1-LRR. The LRR
Molecular docking results of iE-DAP and zNOD1-LRR at different active sites
Grid box
Binding energy (kcal mol1)
Ligand efficiency
No. of H-bonds
Binding coverage area
LRR 1–2 LRR 2–3 LRR 3–4 LRR 4–5 LRR 5–6 LRR 6–7 LRR 7–8 LRR 8–9 N-terminal (a-helix) C-terminal (a-helix) N-terminal (b-sheet) C-terminal (b-sheet) LRR 1–5 (a-helix) LRR 3–7 (a-helix) LRR 5–9 (a-helix) LRR 1–5 (b-sheet) LRR 3–7 (b-sheet) LRR 5–9 (b-sheet) Full a-helix grid Full b-sheet grid Full grid
4.77 3.64 3.57 4.53 3.68 2.46 3.30 2.97 3.51 2.44 4.54 3.21 3.83 3.83 2.39 4.10 3.10 1.07 6.15 4.22 2.72
0.20 0.17 0.16 0.21 0.17 0.11 0.15 0.14 0.16 0.11 0.21 0.15 0.17 0.17 0.11 0.19 0.14 0.05 0.28 0.19 0.12
4 5 4 2 4 4 4 2 8 3 3 4 3 2 5 2 4 2 5 4 3
LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR
1–2 (N-terminal) 2–3 (a-helix) 3–4 (a-helix) 3–5 (a-helix) 4–7 (b-sheet) 8–9 (b-sheet) 8–9 (C-terminal) 9 (C-terminal) 1–2 (a-helix) 7–9 (a-helix) 1–2 (b-sheet) 7–9 (b-sheet) 2–3 (a-helix) 3–4 (a-helix) 6–7 (a-helix) 1–5 (b-sheet) 3–7 (b-sheet) 6–9 (b-sheet) 2–5 (a-helix) 1–5 (b-sheet) 1 (N-terminal)
Bold letters represent ideal binding modes; LRR: leucine rich repeats.
This journal is © The Royal Society of Chemistry 2014
Mol. BioSyst., 2014, 10, 2942--2953 | 2947
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Molecular BioSystems
Paper
regions were defined as the binding sites, represented by different grid box dimensions (see Materials and methods). Out of the 21 candidate docking poses, the pose with iE-DAP bound competently to the a-helix, in the proximity of LRRs 2–5 showed the highest binding score of 6.15 kcal mol1 (Table 2). AutoDock also presented remarkable binding affinities for iE-DAP at other grid sites. To unravel the most probable iE-DAP binding site, six different binding poses were chosen based on their docking score, ligand efficiency, H-bonds, etc. (Table 2). Conceiving the previous NOD1 and iE-DAP interaction reports, these six different complexes (Fig. 5) were considered in a selective fashion to wrap up the complete LRR surfaces and were subjected to MD simulation for stability and conformational analyses. Inference of possible iE-DAP recognition by zNOD1-LRR from MD simulation To elucidate the most probable binding pocket of iE-DAP, six lowest-energy docked conformations (Fig. 5) were simulated for a time period of 30 ns. We performed three independent simulations (10 ns each) for each of the selected docked complex with different random initial velocities in order to acquire an exact ligand–receptor binding mode. The stability of each complex was ascertained by calculating RMSD as a function of simulation time, showing that complex III was the most stable complex as compared to the rest five (Fig. 6). A close analysis of the binding pocket revealed that iE-DAP experienced a noticeable displacement from its starting position in different complexes during MD simulation (Fig. S5, ESI†). In complex I, iE-DAP interacts with the N-terminal loop which is highly flexible in the absence of adjacent NOD1 domains. In complexes II and V, iE-DAP is docked at the convex surface that includes a-helices (a4–8) and connecting loops. During simulation, the
Fig. 6 Root mean square deviation (RMSD) analysis of backbone atoms from 18 independent simulations of 6 different complexes (A) complex I, (B) complex II, (C) complex III, (D) complex IV, (E) complex V and (F) complex VI.
active site is shifted leaving the a-helical binding regions. Complex III displays a well conserved binding pocket where iE-DAP binds to b-sheet (b4–7) and turn regions, which is in agreement with previous mutagenesis data.11 The binding pocket of iE-DAP is well conserved during the simulation. The initial iE-DAP conformation which lies between a7 and a8 in complex IV is dislodged towards the extreme C-termini loop regions.
Fig. 5 Binding site analysis of iE-DAP with zNOD1-LRR using AutoDock 4.2. Protein is shown as cartoon, iE-DAP in stick, and the polar contacts as dotted lines. (A) N-terminal (LRR 1–2); complex I, (B) a-helix (LRR 2–5); complex II, (C) b-sheet (LRR 4–7); complex III, (D) C-terminal (LRR 8–9); complex IV, (E) a-helix (LRR 3–5); complex V and (F) b-sheet (LRR 2–5); complex VI.
2948 | Mol. BioSyst., 2014, 10, 2942--2953
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
Molecular BioSystems
In complex VI, the ligand interacts with the b-sheets at the concave surface (BLRR 2–5). However, the three nitrogen atoms (N8, N9 and N10) in iE-DAP, that were binding to the b-sheet residues before MD, lost the interaction due to a sharp 901 rotation of the iE-DAP from the center of the axis. This testified that complex III represents the most suitable conformation for iE-DAP binding, hence this complex was considered for further biophysical and biochemical analysis. In addition, H-bond analysis of each independent trajectory of complex III revealed that only trajectory I showed constant H-bond interaction with Trp861 (Table 3). Hence, trajectory II and III were discarded by retaining the trajectory I for further evaluation. Apart from H-bond interaction, a close analysis of representative structures (trajectory I of complex III) from each 1 ns showed that iE-DAP constantly formed electrostatic, polar and van der Waals contacts with zNOD1-LRR as summarized in Table 4. For instances, Lys777 and Asp803 were constantly involved in electrostatic interaction whereas, Thr775, Asn831 and Ser833 exhibited polar contacts, and Ile804, Gly805, Trp807, Ser833 and Ile859 formed van der Waals contact during the MD simulation (Table 4). Among the simulated complexes analyzed, the bound iE-DAP obtained a variety of assorted movements and orientations. In each simulation
Table 3 Hydrogen bond interaction in complex III for three different trajectories
Time (ns)
Trajectory I
Trajectory II
Trajectory III
0
Asp803, Trp861 Gly805, Trp861
Asp803, Gly805, Trp861 Thr775, Gly805, Leu832 Trp861 Trp861
Asp803, Gly805, Trp861 Gly805, Trp861
1 2 3
Trp861 Trp807, Trp861 Trp861
4 5
Trp807, Trp861 Trp861 Ser833, Trp861 Ser833, Trp861 Trp861 Asp803, Trp861
6 7 8 9 10
Table 4
Trp807, Trp861 Trp861
— Trp861 Ser833, Trp861
Gly805, Trp861 Trp861
Trp807, Ser833, Trp861 Trp861 Ser833
Ser833, Trp861
Trp861
Ser833, Trp861 —
— Ans831, Ser833, Trp861
Fig. 7 Analysis of intermolecular hydrogen bonds (H-bonds) between iE-DAP and zNOD1-LRR in complex III during 10 ns MD simulation (A) trajectory I, (B) trajectory II and (C) trajectory III.
system, the ligand rotation caused breakage of some crucial H-bonds and other non-bond contacts with the receptor (Fig. S5, ESI†). However, H-bond analyses of complex III revealed that the ligand orientation and displacements were highly restricted, and found to be conserved in all the three trajectories (Fig. 7). A comparative overview of intermolecular interactions between iE-DAP and zNOD1-LRR before and after MD simulations has been summarized in Table S3 (ESI†). Initially, iE-DAP interacted with Thr775, Lys777, Asp803 and Trp861 of zNOD1-LRR through strong H-bonds, indicating the important role of surface exposed polar residues in ligand binding. Interestingly, these residues were found to interact consistently with iE-DAP during 10 ns MD simulations. The visible fluctuations associated with the H-bonds in complex III may be due to the amino acid side chain rotations and different ligand conformations (Fig. 8). The average number of H-bonds calculated as a function of time from these three different trajectories showed approximately equal numbers with no significant variations. Although it is difficult to discriminate
Non-bonded interaction between zNOD1-LRR and iE-DAP at different time scales during 10 ns MD simulation in trajectory I of complex III
Time (ns)
Hydrophobic interaction
Electrostatic interaction
Polar interaction
0 1 2 3 4 5 6 7 8 9 10
Val776 — Trp807 — Trp807, Ile859 — — — — — —
Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777
Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775,
Arg774 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803
Asn831, Asn831 Asn831 Asn831 Asn831, Asn831, Asn831, Asn831 Asn831 Asn831 Asn831
van der Waals interaction Ser833
Ser833 Ser833 Ser833
Ile804, Gly805, Trp807, Trp807, Ile859 Gly805, Ile859 Gly805, Ser833, Ile859 Gly805, Trp807, Ile863 Gly805, Ile859, Ile863 Gly805, Trp807, Ile859 Gly805, Trp807, Ile859, Gly805, Trp807, Ile859, Ile804, Gly805, Trp807, Ile804, Gly805, Trp807,
Ile859
Ile863 Ile863 Ser833, Ile859 Ser833, Ile859
The consistent binding residue Trp861 is shown in bold. This journal is © The Royal Society of Chemistry 2014
Mol. BioSyst., 2014, 10, 2942--2953 | 2949
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Molecular BioSystems
Paper
Fig. 8 Molecular interactions of zNOD1-LRR and iE-DAP complex III in three independent simulations: (A) trajectory I, (B) trajectory II and (C) trajectory III after 10 ns MD simulation.
between the three trajectories based on the number of H-bonds, trajectory I revealed a highly similar H-bonding pattern in comparison to the initial docking conformation. Hence, trajectory I was selected as the best model to describe the ligand recognition by zNOD1-LRR. Fig. 9A and B shows the 3D and 2D
interaction diagrams of zNOD1-LRR and iE-DAP in trajectory I of complex III, respectively. The binding pocket of iE-DAP is comprised of an equal (mixed proportion) distribution of polar and hydrophobic residues and provided an ideal cavity to dock the ligand. A close observation revealed that Trp861 and Trp807
Fig. 9 Intermolecular interactions between iE-DAP and zNOD1-LRR in complex III. (A) Interaction of iE-DAP with the active site amino acids of zNOD1LRR shown in PyMOL. The protein, ligand and hydrogen bonds are shown as cartoon, stick and dotted lines, respectively, (B) 2-dimensional representation of iE-DAP interaction generated using the LigPlot+ program. The protein and ligand are presented in grey and green, respectively. H-bond forming amino acids are shown in grey and hydrophobic contacts as semicircles, (C) the alignment shows the key amino acid residues responsible for iE-DAP and NOD1-LLR interaction in human (red), rohu (green) and zebrafish (blue).
2950 | Mol. BioSyst., 2014, 10, 2942--2953
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
Molecular BioSystems
along with Ile863, Leu832 and Ile859 created a strong hydrophobic environment to anchor the ligand whereas polar residues such as Asp803 and Ser833 hold the ligand through high-affinity H-bond network. The other binding site residues, including Thr775, Lys777, Asn831 and Ser835 contribute crucial electrostatic interactions to the ligand. However, the 2D diagram obtained from LigPlot+67 indicates that the side chain amino group of Asn831 and the side chain protonated nitrogen atom of Trp861 interacted with iE-DAP through H-bonds (Fig. 9B). Overall, it is most likely that the loss of key H-bonds between iE-DAP and zNOD1-LRR sometime during the MD simulation is compensated by electrostatic interactions. Binding free energy calculations The energetic stability of the three trajectories of complex III was further validated by calculating the binding free energies of zNOD1-LRR and iE-DAP complexes using the MM/PBSA method. The binding of zNOD1 and iE-DAP is mostly favored by the coulombic term and is approximately equal in all trajectories of complex III. However, the second and third most contributors (nonpolar term and van der Waals energy) differed significantly. There was an approximation of 1/4th decrease in these contributions in trajectory II and III leading to comparatively stable conformations. Altogether it was found that trajectory I possessed the highest affinity between the binding partners. Therefore, trajectory I was further explored to identify the important residues responsible for the ligand binding. Computational alanine scanning revealed that most of the polar residues interacted with iE-DAP (Thr775, Lys777, Asp803 and Ser833) are crucial for high affinity ligand binding (Table 5). Among the hydrophobic residues, Ile859 was found to be essential for iE-DAP. Trp807 and Trp861, which were engaged in ligand interaction during MD simulation, presented a very negligible effect on the binding free energy. Taken together, it could be stated that zNOD1-LRR mainly recruits polar and negatively charged residues for iE-DAP recognition.
An earlier study showed that the central C-terminal region of the hNOD1 and hNOD2-LRR appeared to be a ‘hot spot’ in terms of muropeptide sensing and this work was well supported by domain swapping and site directed mutagenesis experiments.11 Interestingly, in our previous study it was proposed that LRRs 4–8 of zNOD2-LRR, located centrally to the C-terminal end is responsible for MDP recognition.40 In addition, Girardin et al. in their study showed that His778, Lys790, Gly792, Glu816, Gly818, Trp820 and Trp874 at LRRs 5–8 are indispensable for iE-DAP recognition in hNOD1.11 Similarly a study by Boyle et al. had revealed that the conserved residues on the concave surface of LRR region of hNOD1 and hNOD2 follow a common mechanism of ligand recognition.68 To map the key residues involved in zNOD1–iE-DAP interaction in our study, NOD1 sequences from zebrafish, rohu and human were aligned. As evident from Fig. 9C, the key residues which aid in ligand recognition in both human and zebrafish were conserved and differed in rohu. The reason behind the difference in iE-DAP interacting residues in rohu is not clear and warrants further investigation. In zNOD1, the conserved iE-DAP interacting residues are solvent exposed and located at the concave face of zNOD1-LRR (LRRs 4–7). Among the iE-DAP interacting residues, Lys777, Gly805, Trp807 and Trp861 were conserved in human and zebrafish, while other residues varied. Apart from the conserved residues, three additional residues Asn831, Ser833 and Ile859 were also involved in ligand recognition in zNOD1. In addition, few exceptions were also noticed, for instance, the conserved Gly779 in zebrafish corresponding to Gly792 in human was not found to interact with iE-DAP. Furthermore, Asn831, Ser833 and Ile859 of zebrafish interacted with iE-DAP, but the corresponding residues in human were excluded from the iE-DAP binding pocket. Altogether, our computational approach explained the molecular determinants responsible for ligand recognition in zNOD1. These findings suggested that the conserved LRR domains are highly specific to particular PAMPs in
Table 5 Binding free energy (kJ mol1) calculation of zNOD1-LLR and iE-DAP complexes using the MM/PBSA method: (a) wild zNOD1-LRR complex and (b) mutant zNOD1-LRR complex
Polar contribution
Nonpolar contribution
DGbind
DGcoul
DGps
DGpolar
DGvdW
DGnps
DGnonpolar
(a) Trajectory I Trajectory II Trajectory III
668.00 (5.54) 482.84 (6.03) 364.01 (5.83)
1075.75 (8.54) 1030.09 (7.42) 866.14 (7.41)
618.32 (5.94) 614.88 (5.91) 578.02 (4.72)
457.43 415.21 288.12
197.91 (3.64) 55.53 (4.73) 62.78 (4.43)
12.67 (0.04) 12.10 (0.04) 13.09 (0.03)
210.58 67.63 75.89
(b) Thr775-Ala775 Lys777-Ala777 Asp803-Ala803 Trp807-Ala807 Asn831-Ala831 Ser833-Ala833 Ile859-Ala859 Trp861-Ala861
662.54 665.99 666.82 683.52 675.37 666.45 663.73 680.34
1048.80 (8.55) 1005.69 (8.73) 954.26 (8.26) 1071.53 (8.49) 1079.76 (8.52) 1048.47 (8.49) 1074.79 (8.54) 1048.60 (8.62)
595.20 550.08 498.81 595.54 627.90 593.93 616.96 574.33
453.60 455.61 455.46 475.99 451.86 454.54 457.83 474.28
196.29 197.92 198.87 195.57 198.12 199.16 193.42 206.07
12.66 12.47 12.50 11.95 12.70 12.74 12.49 11.72
208.95 210.39 211.37 207.53 210.82 211.90 205.90 206.07
(5.45) (5.52) (5.34) (5.40) (5.47) (4.43) (5.53) (7.17)
(5.97) (5.88) (5.61) (5.91) (5.91) (5.89) (5.94) (6.10)
(3.71) (3.78) (3.63) (3.45) (3.66) (3.73) (3.71) (3.42)
(0.05) (0.05) (0.04) (0.04) (0.04) (0.05) (0.04) (0.04)
DGbind = binding free energy, DGcoul = electrostatic energy (the coulombic term), DGps = polar solvation energy, DGpolar = polar term (DGcoul + DGps), DGvdW = van der Waals energy, DGnps = nonpolar solvation energy, DGnonpolar = nonpolar term (DGvdW + DGnps), numbers in parentheses indicate standard errors. Experimental binding energy is not available.
This journal is © The Royal Society of Chemistry 2014
Mol. BioSyst., 2014, 10, 2942--2953 | 2951
View Article Online
Molecular BioSystems
different organisms and the ligand binding residues may vary from one species to another.
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Conclusion This study illustrated the modeling of zNOD1-LRR using advanced protein modeling techniques followed by long-range MD simulations. Molecular docking of iE-DAP and zNOD1-LRR followed by MD simulations of the complex indicated that the C-terminal b-sheet region (LRRs 4–7) was responsible for iE-DAP recognition, which agrees with experimental evidences in human. A combination of charged and hydrophobic residues created a compatible binding pocket on zNOD1-LRR for iE-DAP recognition. The residues His775, Lys777, Asp803, Gly805, Trp807, Asn831, Ser833, Ile859 and Trp861 were found to form the binding pocket for iE-DAP recognition. Most of the iE-DAP active site residues in NOD1 were conserved among human and zebrafish, and suggested a strong evolutionary lineage starting from lower to higher eukaryotes. The binding free energy estimation through the MM/PBSA method highlighted the importance of residues His775, Lys777, Asp803, Ser833 and Ile859 for an effective and energetically stable iE-DAP recognition. Consequently, this study unraveled mechanism of iE-DAP recognition in zebrafish, thus providing a benchmark for further investigation on NOD1-based signaling in other vertebrates with a long-term objective of finding novel therapeutic agents to conquer NOD1 related disorders.
Acknowledgements The authors thank the Director, Central Inland Fisheries Research Institute, Barrackpore, Kolkata, India, for providing institutional facility. We also thank Mr Sukanta Kumar Pradhan, Head, Department of Bioinformatics, OUAT, Bhubaneswar, and Dr Sachinandan De, Principal Scientist, NDRI, Karnal, for rendering computational facility and technical advice for improving the manuscript. JM thanks Mr Asim Kumar Jana, Senior Technical Assistant, CIFRI, Barrackpore, Kolkata, India.
References 1 R. Medzhitov and C. A. Janeway, Jr., Semin. Immunol., 1998, 10, 351–353. 2 S. Akira, S. Uematsu and O. Takeuchi, Cell, 2006, 124, 783–801. 3 E. V. Koonin and L. Aravind, Trends Biochem. Sci., 2000, 25, 223–224. ˜ ez, Nat. Rev. Immunol., 2003, 3, 4 N. Inohara and G. Nun 371–382. 5 N. Matsushima, T. Tanaka, P. Enkhbayar, T. Mikami, M. Taga, K. Yamada and Y. Kuroki, BMC Genomics, 2007, 8, 124. 6 K. Hofmann, P. Bucher and J. Tschopp, Trends Biochem. Sci., 1997, 22, 155–156.
2952 | Mol. BioSyst., 2014, 10, 2942--2953
Paper
7 A. M. Verhagen, E. J. Coulson and D. L. Vaux, Genome Biol., 2001, 2, 1–10. 8 J. Bertin and P. S. DiStefano, Cell Death Differ., 2000, 7, 1273–1274. 9 F. Martinon, K. Hofmann and J. Tschopp, Curr. Biol., 2001, 11, 118–120. 10 J. P. Ting, S. B. Willingham and D. T. Bergstralh, Nat. Rev. Immunol., 2008, 8, 372–379. ´hanno, D. Mengin-Lecreulx, P. J. Sansonetti, 11 S. E. Girardin, M. Je P. M. Alzari and D. J. Philpott, J. Biol. Chem., 2005, 280, 38648–38656. 12 M. Chamaillard, S. E. Girardin, J. Viala and D. J. Philpott, Cell. Microbiol., 2003, 5, 581–592. 13 H. Laroui, Y. Yan, Y. Narui, S. A. Ingersoll, S. Ayyadurai, M. A. Charania, F. Zhou, B. Wang, K. Salaita, S. V. Sitaraman and D. Merlin, J. Biol. Chem., 2011, 286, 31003–31013. 14 H. H. Park, Y. C. Lo, S. C. Lin, L. Wang, J. K. Yang and H. Wu, Annu. Rev. Immunol., 2007, 25, 561–586. 15 M. Tohno, W. Ueda, Y. Azuma, T. Shimazu, S. Katoh, J. M. Wang, H. Aso, H. Takada, Y. Kawai, T. Saito and H. Kitazawa, Mol. Immunol., 2008, 45, 194–203. 16 S. H. Oehlers, M. V. Flores, C. J. Hall, S. Swift, K. E. Crosier and P. S. Crosier, Dis. Models & Mech., 2011, 4, 832–841. 17 Z. Sha, J. W. Abernathy, S. Wang, P. Li, H. Kucuktas, H. Liu, E. Peatman and Z. Liu, Dev. Comp. Immunol., 2009, 33, 991–999. 18 W. Q. Chen, Q. Q. Xu, M. X. Chang, P. Nie and K. M. Peng, Fish Shellfish Immunol., 2010, 28, 18–29. 19 B. Swain, M. Basu and M. Samanta, Fish Shellfish Immunol., 2012, 32, 899–908. 20 B. R. Sahoo, B. Swain, M. R. Dikhit, M. Basu, A. Bej, P. Jayasankar and M. Samanta, Appl. Biochem. Biotechnol., 2013, 170, 1282–1309. 21 N. D. Meeker and N. S. Trede, Dev. Comp. Immunol., 2008, 32, 745–757. 22 A. Marchler-Bauer, S. Lu, J. B. Anderson, F. Chitsaz, M. K. Derbyshire, C. DeWeese-Scott, J. H. Fong, L. Y. Geer, R. C. Geer, N. R. Gonzales, M. Gwadz, D. I. Hurwitz, J. D. Jackson, Z. Ke, C. J. Lanczycki, F. Lu, G. H. Marchler, M. Mullokandov, M. V. Omelchenko, C. L. Robertson, J. S. Song, N. Thanki, R. A. Yamashita, D. Zhang, N. Zhang, C. Zheng and S. H. Bryant, Nucleic Acids Res., 2011, 39, D225–D229. 23 E. Quevillon, V. Silventoinen, S. Pillai, N. Harte, N. Mulder, R. Apweiler and R. Lopez, Nucleic Acids Res., 2005, 33, W116–W120. 24 I. Letunic, T. Doerks and P. Bork, Nucleic Acids Res., 2012, 40, D302–D305. 25 R. D. Finn, J. Mistry, J. Tate, P. Coggill, A. Heger, J. E. Pollington, O. L. Gavin, P. Gunasekaran, G. Ceric, K. Forslund, L. Holm, E. L. Sonnhammer, S. R. Eddy and A. Bateman, Nucleic Acids Res., 2010, 38, D211–D222. 26 K. Katoh and D. M. Standley, Mol. Biol. Evol., 2013, 30, 772–780. ¨ffer, R. Agarwala, S. F. Altschul, 27 G. M. Boratyn, A. A. Scha D. J. Lipman and T. L. Madden, Biol. Direct, 2012, 7, 12.
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
28 M. A. Kurowski and J. M. Bujnicki, Nucleic Acids Res., 2003, 31, 3305–3307. 29 O. Fornes, R. Aragues, J. Espadaler, M. A. Marti-Renom, A. Sali and B. Oliva, Bioinformatics, 2009, 25, 1506–1512. 30 Y. Yang, E. Faraggi, H. Zhao and Y. Zhou, Bioinformatics, 2011, 27, 2076–2082. 31 N. Eswar, B. Webb, M. A. Marti-Renom, M. S. Madhusudhan, D. Eramian, M. Y. Shen, U. Pieper and A. Sali, Curr. Protoc. Protein Sci., 2007, ch. 2, Unit 2.9. 32 T. Schwede, J. Kopp, N. Guex and M. C. Peitsch, Nucleic Acids Res., 2003, 31, 3381–3385. 33 C. C. Chen, J. K. Hwang and J. M. Yang, BMC Bioinf., 2009, 10, 366. 34 N. Fernandez-Fuentes, C. J. Madrid-Aliste, B. K. Rai, J. E. Fajardo and A. Fiser, Nucleic Acids Res., 2007, 35, W363–W368. ´onard, X. De Bolle and E. Depiereux, 35 C. Lambert, N. Le Bioinformatics, 2002, 18, 1250–1256. 36 J. Peng and J. Xu, Proteins, 2011, 79, 161–171. 37 S. Wu and Y. Zhang, Nucleic Acids Res., 2007, 35, 3375–3382. 38 Y. Zhang, BMC Bioinf., 2008, 9, 40. ´ll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, 39 S. Pronk, S. Pa M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854. 40 J. Maharana, M. C. Patra, B. C. De, B. R. Sahoo, B. K. Behera, S. De and S. K. Pradhan, J. Mol. Recognit., 2014, 27(5), 260–275. 41 T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092. 42 M. Wiederstein and M. J. Sippl, Nucleic Acids Res., 2007, 35, W407–W410. 43 B. Wallner and A. Elofsson, Protein Sci., 2003, 12, 1073–1086. 44 V. B. Chen, W. B. Arendall, 3rd., J. J. Headd, D. A. Keedy, R. M. Immormino, G. J. Kapral, L. W. Murray, J. S. Richardson and D. C. Richardson, Acta. Crystallogr., Sect. D: Biol. Crystallogr., 2010, 66, 12–21. 45 L. Willard, A. Ranjan, H. Zhang, H. Monzavi, R. F. Boyko, B. D. Sykes and D. S. Wishart, Nucleic Acids Res., 2003, 31, 3316–3319. 46 M. Berjanskii, P. Tang, J. Liang, J. A. Cruz, J. Zhou, Y. Zhou, E. Bassett, C. MacDonell, P. Lu, G. Lin and D. S. Wishart, Nucleic Acids Res., 2009, 37, W670–W677. 47 M. Berjanskii, Y. Liang, J. Zhou, P. Tang, P. Stothard, Y. Zhou, J. Cruz, C. MacDonell, G. Lin, P. Lu and D. S. Wishart, Nucleic Acids Res., 2010, 38, W633–W640. 48 W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
This journal is © The Royal Society of Chemistry 2014
Molecular BioSystems
49 A. K. Malde, L. Zuo, M. Breeze, M. Stroet, D. Poger, P. C. Nair, C. Oostenbrink and A. E. Mark, J. Chem. Theory Comput., 2011, 7, 4026–4037. 50 G. M. Morris, R. Huey, W. Lindstrom, M. F. Sanner, R. K. Belew, D. S. Goodsell and A. J. Olson, J. Comput. Chem., 2009, 30, 2785–2791. 51 D. Spiliotopoulos, A. Spitaleri and G. Musco, PLoS One, 2012, 7, e46902. 52 I. Massova and P. A. Kollman, J. Am. Chem. Soc., 1999, 121, 8133–8143. 53 N. A. Baker, D. Sept, S. Joseph, M. J. Holst and J. A. McCammon, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 10037–10041. 54 S. P. Brown and S. W. Muchmore, J. Med. Chem., 2009, 52, 3159–3165. 55 A. V. Kajava and B. Kobe, Protein Sci., 2002, 11, 1082–1090. 56 J. Maharana, B. Swain, B. R. Sahoo, M. R. Dikhit, M. Basu, A. S. Mahapatra, P. Jayasankar and M. Samanta, Fish Physiol. Biochem., 2012, 39, 1007–1023. 57 B. R. Sahoo, J. Maharana, G. K. Bhoi, S. K. Lenka, M. C. Patra, M. R. Dikhit, P. K. Dubey, S. K. Pradhan and B. K. Behera, Mol. Biosyst., 2014, 10, 1104–1116. 58 D. Baker and A. Sali, Science, 2001, 294, 93–96. 59 B. Kobe and J. Deisenhofer, Nature, 1995, 374, 183–186. 60 R. J. Johnson, J. G. McCoy, C. A. Bingman, G. N. Phillips Jr. and R. T. Raines, J. Mol. Biol., 2007, 368, 434–449. 61 D. W. Buchan, F. Minneci, T. C. Nugent, K. Bryson and D. T. Jones, Nucleic Acids Res., 2013, 41, W349–W357. 62 G. N. Ramachandran, C. Ramakrishnan and V. Sasisekharan, J. Mol. Biol., 1963, 7, 95–99. 63 R. A. Laskowski, M. W. MacArthur, D. S. Moss and J. M. Thornton, J. Appl. Crystallogr., 1993, 26, 283–291. ¨thy, J. U. Bowie and D. Eisenberg, Nature, 1992, 356, 64 R. Lu 83–85. 65 C. Colovos and T. O. Yeates, Protein Sci., 1993, 2, 1511–1519. 66 T. Tanabe, M. Chamaillard, Y. Ogura, L. Zhu, S. Qiu, J. Masumoto, P. Ghosh, A. Moran, M. M. Predergast, ˜ ez, ´n G. Tromp, C. J. Williams, N. Inohara and G. Nu EMBO J., 2004, 23, 1587–1597. 67 R. A. Laskowski and M. B. Swindells, LigPlot+: multiple ligand-protein interaction diagrams for drug discovery, J. Chem. Inf. Model., 2011, 51, 2778–2786. 68 J. P. Boyle, S. Mayle, R. Parkhouse and T. P. Monie, Front. Microb. Immunol., 2013, 4, 317.
Mol. BioSyst., 2014, 10, 2942--2953 | 2953
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
PAPER
Cite this: Mol. BioSyst., 2014, 10, 2942
View Journal | View Issue
Structural and functional investigation of zebrafish (Danio rerio) NOD1 leucine rich repeat domain and its interaction with iE-DAP† Jitendra Maharana,*a Bikash Ranjan Sahoo,b Aritra Bej,b Mahesh Chandra Patra,bc Budheswar Dehury,d Gopal Krushna Bhoi,b Santosh Kumar Lenka,b Jyoti Ranjan Sahoo,b Ajaya Kumar Routb and Bijay Kumar Beheraa Nucleotide binding and oligomerization domain 1 (NOD1), a cytoplasmic pattern recognition receptor (PRR) and is a key component for modulating innate immunity and signaling. It is highly specific to g-D-Glu-mDAP (iE-DAP), a cell wall component of Gram-negative and few Gram-positive bacteria. In the absence of the experimental structure of NOD1 leucine rich repeat (NOD1-LRR) domain, the NOD signaling cascade mediated through NOD1 and iE-DAP interaction is poorly understood. Herein, we modeled 3D structure of zebrafish NOD1-LRR (zNOD1-LRR) through a protein-threading approach and structural integrity of the model was assessed using molecular dynamics simulations. Molecular interaction analysis of iE-DAP and zNOD1-LRR, their complex stability and binding free energy studies were conducted to anticipate the ligand binding residues in zNOD1. Our study revealed that His775, Lys777, Asp803, Gly805, Trp807, Asn831, Ser833, Ile859 and Trp861 situated in the b-sheet region of zNOD1LRR could be involved in iE-DAP recognition, which correlates the earlier findings in human. Comparison
Received 3rd April 2014, Accepted 29th July 2014
of binding free energies of native and mutant zNOD1–iE-DAP complexes delineated His775, Lys777,
DOI: 10.1039/c4mb00212a
This study provides the first comprehensive description of biophysical and biochemical parameters
Asp803, Ser833 and Ile859 as the pivotal residues for energetic stability of NOD1 and iE-DAP interaction. responsible for NOD1 and iE-DAP interaction in zebrafish, which is expected to shed more light on
www.rsc.org/molecularbiosystems
NOD1 signaling and therapeutic applications in other organisms.
Introduction Innate immune system is the most ancestral and ubiquitous host defense mechanism to combat a wide range of pathogens. It is activated upon specific recognition of pathogen-associated molecular patterns (PAMPs) by pattern recognition receptors (PRRs).1,2 PRRs are widespread in extracellular, membrane and cytoplasmic regions, and are classified according to their PAMP specificity, function, localization and evolutionary origin. The nucleotide binding and oligomerization domain (NOD)-like
a
Biotechnology Laboratory, Central Inland Fisheries Research Institute, Barrackpore, Kolkata-700120, West Bengal, India. E-mail: [email protected]; Tel: +91 3325921190/91 b BIF-Centre, Department of Bioinformatics, Orissa University of Agriculture and Technology, Bhubaneswar-751003, Odisha, India c Animal Genomics Laboratory, Animal Biotechnology Centre, National Dairy Research Institute, Karnal-132001, Haryana, India d Department of Life Science and Bioinformatics, Assam University, Silchar-788011, Assam, India † Electronic supplementary information (ESI) available. See DOI: 10.1039/ c4mb00212a
2942 | Mol. BioSyst., 2014, 10, 2942--2953
receptors (NLRs) are some of the various reported PRRs, which belongs to the family of intracellular PRRs having canonical tripartite domain architecture. This includes an N-terminal effector binding domain (EBD), a central NACHT domain mediating self-regulation and oligomerization3 and a C-terminal leucine rich repeat (LRR) domain that recognizes PAMPs.4,5 The EBD domain is composed of either one or two caspase activating and recruitment domains (CARD),6 or a baculovirus inhibitor of apoptosis protein repeat (BIR),7 or a pyrin domain (PYD).8,9 Based on the diversity of the N-terminal EBD, NLRs are divided into five major subgroups, viz., NLRA (NLR family, acidic domain containing proteins), NLRB (NACHT, LRR and BIR-containing proteins), NLRC (NACHT, LRR and CARD-containing proteins), NLRP (NACHT, LRR and PYD-containing proteins) and NLRX (NLR family with no strong homology to the N-terminal domain of any other NLR subfamily member).10 The NLRC or NOD sub-group is constituted of five receptor proteins, NLRC1 (NOD1), NLRC2 (NOD2), NLRC3 (NOD3), NLRC4 (IPAF) and NLRC5 (NOD4), which sense various microbial derivatives and activate a complex downstream signaling cascade. Among these, NOD1 recognizes peptidoglycan (PGN)
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
derivatives, such as g-D-Glu-mDAP (iE-DAP), produced by all Gram-negative and specific Gram-positive bacteria.11–13 The NOD–PAMP interaction involves ATP recruitment which further leads to oligomerization of the NACHT domain, followed by CARD–CARD interaction between NOD and receptor-interacting serine/threonine-protein kinase 2 (RIP2), that engenders NF-kB activation and signal transduction.14 Previous studies showed that, NOD1 plays a significant role in innate immunity of higher eukaryotes.12,13,15 NOD1 has been recently characterized in lower eukaryotes like fish, frogs, snakes, etc. Notably, the role of the Nod1 gene in host protection has been studied in several bony fishes including zebrafish (Danio rerio),16 catfish (Ictalurus punctatus),17 grass carp (Ctenopharyngodon idella)18 and rohu (Labeo rohita).19,20 Over the past decade, with the advent of next-generation sequencing technologies, considerable advancements have been made in unraveling the mechanism of host–pathogen interactions. But, lack of structural information on NLR family members possesses a constraint to understand pathogen recognition and immune signaling pathways at the molecular level in higher-eukaryotes. In this scenario, zebrafish is recognized as an attractive model organism to study the in vivo importance of NLR proteins. In a recent study, Meeker et al.21 have shown that zebrafish shares notable resemblance with the human immune system. Furthermore, a study of Oehlers et al.16 has revealed that NLR protein family members, particularly NOD1 and NOD2 in zebrafish, are responsible for inflammatory bowel disease (IBD) predisposition. In this study, we have constructed the three dimensional (3D) model of zebrafish NOD1-LRR (zNOD1-LRR), using various molecular modeling techniques that includes multiple-template homology modeling, protein threading and ab initio modeling. We have established the low energy binding modes of iE-DAP and elucidated intermolecular interactions between iE-DAP and zNOD1-LRR using molecular docking, followed by molecular dynamics (MD) simulations. The binding free energy between zNOD1-LRR and iE-DAP was calculated using the Molecular Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) method. Taken together, this study for the first time demonstrates structural and dynamic properties of zNOD1-LRR in association with iE-DAP.
Materials and methods Domain analysis The protein sequence of zNOD1 was retrieved from the NCBI protein database (XP_002665106). The Conserved Domain Database (CDD) was used for analyzing the domain architecture of zNOD1.22 In addition to this, the zNOD1 protein sequence was also scanned in InterProScan,23 SMART24 and Pfam25 database to locate N-terminal CARD, central NACHT and C-terminal LRR domains. The LRR regions were manually identified by aligning the zNOD1-LRR domain with human, mouse and rohu. The multiple sequence alignment was executed in MAFFT26 and LRR regions were manually identified by locating ‘‘LxxLxLxxNxL’’
This journal is © The Royal Society of Chemistry 2014
Molecular BioSystems
motif (where L = leucine/isoleucine/valine/phenylalanine; x = any amino acid; N = asparagine/threonine/serine/cysteine). Template identification and molecular modeling To find suitable templates for the modeling of zNOD1-LRR, DELTA-BLAST27 search was performed against the protein data bank (PDB) (http://www.pdb.org/). Because of low sequence identity between target protein and PDB templates, the protein sequence was submitted to various protein threading servers, such as Genesilico Metaserver2,28 ModLink+29 and SPARKS-X,30 to identify suitable templates through fold recognition methods. Three best scoring templates resulting from fold recognition servers were considered for model construction using MODELLER 9.12.31 We have adopted both single-template and multipletemplate approaches for building zNOD1-LRR structures. The 3D models were also constructed using the automated modeling programs, SWISS-MODEL,32 BioSerf, (PS)2,33 M4T,34 EasyPred3D,35 Raptor-X,36 LOMETS37 and I-TASSER38 for comparison. The MODELLER-derived 3D models of zNOD1-LRR were evaluated based on discrete optimized protein energy (DOPE) scores. Three models with lowest DOPE scores were chosen for further refinement using GROMACS 4.5.5 simulation package.39 Molecular dynamics simulations MD simulations of the zNOD1-LRR models were carried out using Gromos96-53a6 force field and SPC/E water models in separate cubic boxes with a minimum distance of 15 Å between the protein surfaces and box edges. To neutralize the simulation systems, a physiological ionic strength (0.15 M NaCl) of counter ions was added. The atomic composition of the simulation systems is provided in Table S1 (ESI†). The final MD simulation of 20 ns was carried out for each system and the simulation parameters were followed from the previous report.40 The particle mesh Ewald (PME) summation method41 was used to treat the long-range electrostatic interactions. The non-bonded interactions were truncated after 12 Å of the cut-off radius. The temperature was coupled to 300 K using the Berendsen thermostat algorithm. A time step of 2 fs was used and the non-bonded neighbor list was updated heuristically. Conformations were saved at 2 ps time intervals for analysis. The average structure was extracted from the most stable portion of dynamics trajectory followed by energy minimization. The best model was chosen based on the stability of backbone atoms, individual residues, radius of gyration (Rg), secondary structure conservation and residual fluctuation. Model validation The stereochemical quality of the final energy-optimized model was verified using SAVeS (http://nihserver.mbi.ucla.edu/SAVES/), ProSA42 and ProQ.43 Bond length and bond angle analyses of the model was performed in MolProbity.44 Furthermore, the Z-score of hydrogen bond (H-bond) energy, packing defects, bump score, radius of gyration and deviation of Y angles of the refined model was tested in VADAR,45 GeNMR46 and PROSESS47 web servers. Time-dependent secondary structure analyses of the models were executed using visual molecular dynamics (VMD 1.9.1).48
Mol. BioSyst., 2014, 10, 2942--2953 | 2943
View Article Online
Molecular BioSystems
Paper
Two dimensional (2D) plots were generated using Grace 5.1.23 programs (http://plasma-gate.weizmann.ac.il/Grace/).
EMM = Eint + Ecoul + EvdW
Molecular docking of iE-DAP
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
The molecular mechanics interaction energy (EMM), is defined as (eqn (3)):
The 2D structure of iE-DAP (CID_45480617) was obtained from the PubChem database (http://pubchem.ncbi.nlm.nih.gov/). The Automated Topology Builder (ATB) server49 was used to generate the 3D coordinates of iE-DAP employing energy minimization with the Gromos96-53a6 force field. Molecular docking was carried out using AutoDock 4.250 to identify key ligand binding residues. The zNOD1-LRR model was assigned with Kollman charges and Geister partial charges were applied for iE-DAP using Autodock Tools 1.5. The binding site information of zNOD1-LRR was obtained from previously published literature.11,20 Initially, blind docking was performed to identify probable iE-DAP interaction sites. The docking was performed with different grid dimension around various LRR regions with a grid spacing of 0.375 Å. The initial population size was 300 and maximum number energy of evaluation was set to 2.5 107. The docked conformations were clustered with a 2 Å cut-off root mean square deviation (RMSD), and the resulting docked poses were ranked according to their binding energy scores. The best complexes from each docking simulation were subjected to MD simulations to optimize the ligand–receptor interaction.
where Eint represents bond, angle and torsion angle energies; Ecoul stands for electrostatic energy; and EvdW denotes van der Waals energy. The solvation free energy term (Gsol), is calculated using eqn (4): Gsol = Gpolar + Gnonpolar
A total of six docked zNOD1-LRR–iE-DAP complexes were selected on the basis of binding energy scores and the number of intermolecular H-bonds. Three independent MD simulations were performed on the selected docking complexes for better sampling of results. The ligand topology was prepared using ATB server.49 The procedure and parameters of MD simulations were the same as described in our previous work.40 Position restraints were applied to the backbone atoms of the modeled protein during equilibration simulation lasting 100 ps with both NVT and subsequent NPT conditions. The temperature was weakly coupled at 300 K via Berendsen method. The pressure was maintained at 1 bar through Parrinello–Rahman barostat algorithm in the NPT ensemble. After temperature and pressure equilibration, a total of 10 ns production MD was carried out for each complex. To amplify the performance and scaling of production recording, each simulation was truncated with three equal time intervals. For all six complexes, a total of 18 truncated MD (3 for each complex) was performed and the trajectories were analyzed. Binding free energy calculation A total of 500 snapshots were obtained from the MD trajectories to calculate the binding free energy51 by employing the MM/PBSA method52 as described in eqn (1): DGbinding = Gcomplex (GzNOD1 + GiE-DAP)
(1)
The individual energy terms for protein (zNOD1), ligand (iE-DAP) and complex are computed based on eqn (2): hGi = hEMMi + hGsoli ThSMMi
2944 | Mol. BioSyst., 2014, 10, 2942--2953
(2)
(4)
where Gpolar and Gnonpolar terms were calculated using APBS program.53 The polar term (Gpolar) was calculated by solving the non-linearized PB equation. The parameters for APBS calculation are as follows: grid spacing, 0.5 Å; temperature, 296 K; and salt concentration, 0.15 M. The nonpolar contribution, Gnonpolar, was computed using eqn (5): Gnonpolar = gSASA + b
(5)
where g = 0.0227 kJ mol1 Å2 and b = 0 kJ mol1.54 The dielectric boundary was defined using a probe radius of 1.4 Å. The standard errors were calculated using the following eqn (6): Standard error (SE) = s/ON
MD simulation of complexes
(3)
(6)
where s is the standard deviation and N is the number of structures used in the calculation.
Results and discussion Domain identification The domain analysis revealed that, zNOD1 sequence (940 amino acids (aa)) is comprised of a short N-terminal CARD (13–103 aa), a central NACHT domain (190–525 aa) and a C-terminal LRR domain (689–933 aa). Multiple sequence alignment of NOD1 sequences of zebrafish, human, mouse and rohu delineating the functional domains is shown in Fig. S1 (ESI†). The LRR domain in zNOD1 consisted of a series of covalently linked 21–29 aa residue motifs with the conserved pattern ‘‘LxxLxLxxNxL’’. Structurally, each LRR motif is comprised of one b-sheet with the pattern ‘xxLxLxx’ and one a-helix connected by a loop.55 Automatic prediction of LRR regions using SMART, Pfam and InterProScan advocated six LRR motifs in zNOD1. However, manual identification of LRR motifs from the aligned sequences of NOD1 from different species showed three additional LRR motifs (Fig. 1).20,40,56,57 Molecular modeling of zNOD1-LRR The search for templates for zNOD1-LRR in PDB database showed a very low percentage of sequence identities (i.e., o30%; twilight zone).58 So, the modeling exercises were performed employing fold recognition or threading methods.30 The resultant templates suggested by different threading programs are shown in Table S2 (ESI†). The 3D coordinates of zNOD1-LRR was assigned using crystal structures of porcine (PDB ID: 2BNH59),
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
Molecular BioSystems
Fig. 1 Multiple sequence alignment of human, mouse, rohu and zebrafish NOD1-LRR sequences using MAFFT. The leucine rich regions (LRRs) of zebrafish NOD1 are highlighted in boxes and the nine LRRs motifs are shown at the top of the alignment. The symbols ‘*’, ‘:’ and ‘.’ represents identical, conserved and semi-conserved substitutions, respectively.
mouse (3TSR) and human (1Z7X60) ribonuclease inhibitors. Models obtained from the automated molecular modeling servers were compared with those of MODELLER. To investigate the accuracy of built models, PSIPRED61 generated secondary structures of zNOD1-LRR were compared with the above mentioned template structures (Fig. S2, ESI†). After careful inspection, three top scored models were selected [model-1 (multiple templates); model-2 (Raptor-X); and model-3 (I-TASSER)] based on their DOPE scores (Fig. S3, ESI†). Analysis of MD trajectories The 20 ns MD trajectories of the selected zNOD1-LRR models were studied to gain insight into the structural characteristics. First, RMSD of backbone atoms was calculated as a function of simulation time. Both model-1 and 2 showed an average backbone deviation of 3.0 Å from its initial conformation; however, model-3 experienced an unusual backbone deviation (B6.5 Å) between 10 and 20 ns (Fig. 2A). The overall stability of zNOD1LRR was further confirmed by calculating Rg as a function of simulation time. The constant gyradius of 19.5 Å suggested a consistent shape and size of model-1 during simulation (Fig. 2B). In addition, the root mean square fluctuations (RMSF) of Ca atoms were calculated for each residue. The Ca atoms of
This journal is © The Royal Society of Chemistry 2014
Fig. 2 Evaluation of stability parameters of the modeled zNOD1-LRR in the three independent MD simulations (20 ns). (A) Root mean square deviations (RMSD), (B) radius of gyration (Rg) and (C) root mean square fluctuations (RMSFs). The color legends for different models are shown above the plot.
Mol. BioSyst., 2014, 10, 2942--2953 | 2945
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Molecular BioSystems
terminal amino acids fluctuated up to B10 Å in all models. The ligand binding domain C-terminal as studied earlier showed a prominent motion in model-1. LRR 6–8 regions in model-1 presented a distinct flipped conformation and indicated its engagement in ligand binding. On the other hand in model-2 and 3 these regions were comparatively stable (Fig. 2C). The C-terminal loop and helical regions of LRR domain were flexible during 20 ns simulation and suggested a structural rearrangement of functional domains in NOD1 upon PAMPs recognition.20 Secondary structure evolution from the trajectories showed that in model-1 all nine a-helices and nine b-strands retained their secondary structural properties throughout the MD simulation in aqueous solution (Fig. 3). Moreover, conformational dynamics analysis revealed no significant conformational alterations in model-1 during MD simulations as compared to model-2 and 3 (Fig. S4, ESI†). All together the MD simulation of top three models advocated model-1 as the most stable conformation and hence considered for further structural analysis. Structural analysis of the zNOD1-LRR model The optimized model obtained after MD simulation was investigated to delineate its key LRR domain and spatial conformation. The modeled zNOD1-LRR was displayed as a semicircle with a concave surface containing nine parallel b-strands, a convex face constituting nine a helices and short loops connecting these a/b-subunits. The LRRs within zNOD1-LRR are well conserved across different species as revealed from the sequence
Paper
Fig. 4 Schematic representation of the zNOD1-LRR model showing different LRR regions. The a-helices are shown in red, b-sheets in yellow and turns in green. Alignment of nine LRR motifs with the conserved pattern ‘‘LxxLxLxxNxL’’ is shown. The a-helices and b-sheets are displayed in red and blue charterers, respectively. The iE-DAP binding sites evidenced from previous reports are shown at the right bottom (transparent surface representations).
alignment (Fig. 4). Our study demonstrated that, the ligand binding active site of zNOD1 lies within LRRs 4–7 of the b-sheet region and corresponded to the previous reports in human (LRRs 5–7) by Girardin et al.11
Fig. 3 Secondary structures of three different models (A) model-1, (B) model-2 and (C) model-3 as a function of simulation time. Color segments given at the bottom of the figure represents the different secondary structural properties.
2946 | Mol. BioSyst., 2014, 10, 2942--2953
This journal is © The Royal Society of Chemistry 2014
View Article Online
Paper
Molecular BioSystems
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Model validation The stereo-chemical properties of the proposed zNOD1-LRR model were validated using various protein structure validation servers. The accuracy of dihedral angles (F/C) was gauzed by drawing the Ramachandran plot62 in PROCHECK integrated at SAVeS server. The zNOD1-LRR model showed that 77.90 and 98.60% of residues fell in the most favored and allowed regions, respectively, with an overall G-factor score of 0.04 (cut-off score 40.50) (Table 1). However, 1.40% of residues that includes Lys743, Gln864 and Ser868 fell within the disallowed regions. A close inspection of the modeled structures and previous reports revealed that the disallowed residues neither Table 1
Model validation report of zNOD1-LRR
Servers
zNOD1-LRR
PROCHECK
Verify3D ERRAT ProSA ProQ MolProbity VADAR
GeNMR PROSESS
Table 2
Most favored regions (%) Additionally allowed regions (%) Generously allowed regions (%) Disallowed regions (%) Overall G-factor Averaged 3D–1D score 40.2 Overall quality Z-Score LGscore MaxSub Cb deviations 40.25 Å Residues with bad bonds (%) Residues with bad angles (%) Standard deviation of w1 pooled Mean H-bond energy Generously allowed O angles (%) Packing defects (%) Percentage of 95% buried residues Ramachandran outside of most favored Bump score Radius gyration score Deviation of Y angles w1 score
77.90 18.90 1.80 1.40 0.04 99.59 75.10 7.09 6.33 0.54 0.00 0.00 0.00 2.77 0.70 0.26 0.21 0.54 2.56 0.34 1.20 1.23 1.20
involved in ligand binding nor constituted the structurally rigid regions. This indicated that the stereo-chemical properties of the model were within the acceptable range.63 The Verify3D score of the proposed model was 99.59%, which is above the cut off (Z80%) value and indicated that the quality of the proposed model is reasonably good.64 The ERRAT score provides accuracy of the non-bonded atomic contacts and our predicted model has a score of 75.10%, which is greater than the acceptable value (50%).65 The Z-score of our proposed model calculated using ProSA was in agreement with that of the experimental PDB structures of similar sizes. ProQ analysis indicated that the quality of zNOD1-LRR model was ‘extremely good’ based on their MaxSub and LGscores. Bond length and bond angle analysis of the proposed model using MolProbity revealed that none of the residues had bad side chain or main chain conflicts. The H-bond energy, packing defect, bump score, Rg and deviation of Y angles of the proposed model as analyzed by VADAR, GeNMR and PROSESS servers were within the cutoff range. A detailed validation report is summarized in Table 1. Analysis of docking results To identify the possible iE-DAP binding site in zNOD1-LRR, several molecular docking simulations were performed. Earlier studies have revealed that in NOD1 and NOD2, the b-sheet regions of the C-terminus (hot spot) appears to be critical for muropeptide recognition.11,40,56,66 In addition, a recent computational study on rohu NOD1 predicted C-terminal LRRs 8–9 as the iE-DAP binding site.20 Keeping all the evidences in support of the importance of C-terminal conserved LRR domain for iE-DAP recognition encouraged us to further investigate the binding modes between zNOD1-LRR and iE-DAP. In this study, we performed a total of 21 independent docking calculations to investigate the interaction of iE-DAP with zNOD1-LRR. The LRR
Molecular docking results of iE-DAP and zNOD1-LRR at different active sites
Grid box
Binding energy (kcal mol1)
Ligand efficiency
No. of H-bonds
Binding coverage area
LRR 1–2 LRR 2–3 LRR 3–4 LRR 4–5 LRR 5–6 LRR 6–7 LRR 7–8 LRR 8–9 N-terminal (a-helix) C-terminal (a-helix) N-terminal (b-sheet) C-terminal (b-sheet) LRR 1–5 (a-helix) LRR 3–7 (a-helix) LRR 5–9 (a-helix) LRR 1–5 (b-sheet) LRR 3–7 (b-sheet) LRR 5–9 (b-sheet) Full a-helix grid Full b-sheet grid Full grid
4.77 3.64 3.57 4.53 3.68 2.46 3.30 2.97 3.51 2.44 4.54 3.21 3.83 3.83 2.39 4.10 3.10 1.07 6.15 4.22 2.72
0.20 0.17 0.16 0.21 0.17 0.11 0.15 0.14 0.16 0.11 0.21 0.15 0.17 0.17 0.11 0.19 0.14 0.05 0.28 0.19 0.12
4 5 4 2 4 4 4 2 8 3 3 4 3 2 5 2 4 2 5 4 3
LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR LRR
1–2 (N-terminal) 2–3 (a-helix) 3–4 (a-helix) 3–5 (a-helix) 4–7 (b-sheet) 8–9 (b-sheet) 8–9 (C-terminal) 9 (C-terminal) 1–2 (a-helix) 7–9 (a-helix) 1–2 (b-sheet) 7–9 (b-sheet) 2–3 (a-helix) 3–4 (a-helix) 6–7 (a-helix) 1–5 (b-sheet) 3–7 (b-sheet) 6–9 (b-sheet) 2–5 (a-helix) 1–5 (b-sheet) 1 (N-terminal)
Bold letters represent ideal binding modes; LRR: leucine rich repeats.
This journal is © The Royal Society of Chemistry 2014
Mol. BioSyst., 2014, 10, 2942--2953 | 2947
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Molecular BioSystems
Paper
regions were defined as the binding sites, represented by different grid box dimensions (see Materials and methods). Out of the 21 candidate docking poses, the pose with iE-DAP bound competently to the a-helix, in the proximity of LRRs 2–5 showed the highest binding score of 6.15 kcal mol1 (Table 2). AutoDock also presented remarkable binding affinities for iE-DAP at other grid sites. To unravel the most probable iE-DAP binding site, six different binding poses were chosen based on their docking score, ligand efficiency, H-bonds, etc. (Table 2). Conceiving the previous NOD1 and iE-DAP interaction reports, these six different complexes (Fig. 5) were considered in a selective fashion to wrap up the complete LRR surfaces and were subjected to MD simulation for stability and conformational analyses. Inference of possible iE-DAP recognition by zNOD1-LRR from MD simulation To elucidate the most probable binding pocket of iE-DAP, six lowest-energy docked conformations (Fig. 5) were simulated for a time period of 30 ns. We performed three independent simulations (10 ns each) for each of the selected docked complex with different random initial velocities in order to acquire an exact ligand–receptor binding mode. The stability of each complex was ascertained by calculating RMSD as a function of simulation time, showing that complex III was the most stable complex as compared to the rest five (Fig. 6). A close analysis of the binding pocket revealed that iE-DAP experienced a noticeable displacement from its starting position in different complexes during MD simulation (Fig. S5, ESI†). In complex I, iE-DAP interacts with the N-terminal loop which is highly flexible in the absence of adjacent NOD1 domains. In complexes II and V, iE-DAP is docked at the convex surface that includes a-helices (a4–8) and connecting loops. During simulation, the
Fig. 6 Root mean square deviation (RMSD) analysis of backbone atoms from 18 independent simulations of 6 different complexes (A) complex I, (B) complex II, (C) complex III, (D) complex IV, (E) complex V and (F) complex VI.
active site is shifted leaving the a-helical binding regions. Complex III displays a well conserved binding pocket where iE-DAP binds to b-sheet (b4–7) and turn regions, which is in agreement with previous mutagenesis data.11 The binding pocket of iE-DAP is well conserved during the simulation. The initial iE-DAP conformation which lies between a7 and a8 in complex IV is dislodged towards the extreme C-termini loop regions.
Fig. 5 Binding site analysis of iE-DAP with zNOD1-LRR using AutoDock 4.2. Protein is shown as cartoon, iE-DAP in stick, and the polar contacts as dotted lines. (A) N-terminal (LRR 1–2); complex I, (B) a-helix (LRR 2–5); complex II, (C) b-sheet (LRR 4–7); complex III, (D) C-terminal (LRR 8–9); complex IV, (E) a-helix (LRR 3–5); complex V and (F) b-sheet (LRR 2–5); complex VI.
2948 | Mol. BioSyst., 2014, 10, 2942--2953
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
Molecular BioSystems
In complex VI, the ligand interacts with the b-sheets at the concave surface (BLRR 2–5). However, the three nitrogen atoms (N8, N9 and N10) in iE-DAP, that were binding to the b-sheet residues before MD, lost the interaction due to a sharp 901 rotation of the iE-DAP from the center of the axis. This testified that complex III represents the most suitable conformation for iE-DAP binding, hence this complex was considered for further biophysical and biochemical analysis. In addition, H-bond analysis of each independent trajectory of complex III revealed that only trajectory I showed constant H-bond interaction with Trp861 (Table 3). Hence, trajectory II and III were discarded by retaining the trajectory I for further evaluation. Apart from H-bond interaction, a close analysis of representative structures (trajectory I of complex III) from each 1 ns showed that iE-DAP constantly formed electrostatic, polar and van der Waals contacts with zNOD1-LRR as summarized in Table 4. For instances, Lys777 and Asp803 were constantly involved in electrostatic interaction whereas, Thr775, Asn831 and Ser833 exhibited polar contacts, and Ile804, Gly805, Trp807, Ser833 and Ile859 formed van der Waals contact during the MD simulation (Table 4). Among the simulated complexes analyzed, the bound iE-DAP obtained a variety of assorted movements and orientations. In each simulation
Table 3 Hydrogen bond interaction in complex III for three different trajectories
Time (ns)
Trajectory I
Trajectory II
Trajectory III
0
Asp803, Trp861 Gly805, Trp861
Asp803, Gly805, Trp861 Thr775, Gly805, Leu832 Trp861 Trp861
Asp803, Gly805, Trp861 Gly805, Trp861
1 2 3
Trp861 Trp807, Trp861 Trp861
4 5
Trp807, Trp861 Trp861 Ser833, Trp861 Ser833, Trp861 Trp861 Asp803, Trp861
6 7 8 9 10
Table 4
Trp807, Trp861 Trp861
— Trp861 Ser833, Trp861
Gly805, Trp861 Trp861
Trp807, Ser833, Trp861 Trp861 Ser833
Ser833, Trp861
Trp861
Ser833, Trp861 —
— Ans831, Ser833, Trp861
Fig. 7 Analysis of intermolecular hydrogen bonds (H-bonds) between iE-DAP and zNOD1-LRR in complex III during 10 ns MD simulation (A) trajectory I, (B) trajectory II and (C) trajectory III.
system, the ligand rotation caused breakage of some crucial H-bonds and other non-bond contacts with the receptor (Fig. S5, ESI†). However, H-bond analyses of complex III revealed that the ligand orientation and displacements were highly restricted, and found to be conserved in all the three trajectories (Fig. 7). A comparative overview of intermolecular interactions between iE-DAP and zNOD1-LRR before and after MD simulations has been summarized in Table S3 (ESI†). Initially, iE-DAP interacted with Thr775, Lys777, Asp803 and Trp861 of zNOD1-LRR through strong H-bonds, indicating the important role of surface exposed polar residues in ligand binding. Interestingly, these residues were found to interact consistently with iE-DAP during 10 ns MD simulations. The visible fluctuations associated with the H-bonds in complex III may be due to the amino acid side chain rotations and different ligand conformations (Fig. 8). The average number of H-bonds calculated as a function of time from these three different trajectories showed approximately equal numbers with no significant variations. Although it is difficult to discriminate
Non-bonded interaction between zNOD1-LRR and iE-DAP at different time scales during 10 ns MD simulation in trajectory I of complex III
Time (ns)
Hydrophobic interaction
Electrostatic interaction
Polar interaction
0 1 2 3 4 5 6 7 8 9 10
Val776 — Trp807 — Trp807, Ile859 — — — — — —
Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777, Lys777
Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775, Thr775,
Arg774 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803 Asp803
Asn831, Asn831 Asn831 Asn831 Asn831, Asn831, Asn831, Asn831 Asn831 Asn831 Asn831
van der Waals interaction Ser833
Ser833 Ser833 Ser833
Ile804, Gly805, Trp807, Trp807, Ile859 Gly805, Ile859 Gly805, Ser833, Ile859 Gly805, Trp807, Ile863 Gly805, Ile859, Ile863 Gly805, Trp807, Ile859 Gly805, Trp807, Ile859, Gly805, Trp807, Ile859, Ile804, Gly805, Trp807, Ile804, Gly805, Trp807,
Ile859
Ile863 Ile863 Ser833, Ile859 Ser833, Ile859
The consistent binding residue Trp861 is shown in bold. This journal is © The Royal Society of Chemistry 2014
Mol. BioSyst., 2014, 10, 2942--2953 | 2949
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Molecular BioSystems
Paper
Fig. 8 Molecular interactions of zNOD1-LRR and iE-DAP complex III in three independent simulations: (A) trajectory I, (B) trajectory II and (C) trajectory III after 10 ns MD simulation.
between the three trajectories based on the number of H-bonds, trajectory I revealed a highly similar H-bonding pattern in comparison to the initial docking conformation. Hence, trajectory I was selected as the best model to describe the ligand recognition by zNOD1-LRR. Fig. 9A and B shows the 3D and 2D
interaction diagrams of zNOD1-LRR and iE-DAP in trajectory I of complex III, respectively. The binding pocket of iE-DAP is comprised of an equal (mixed proportion) distribution of polar and hydrophobic residues and provided an ideal cavity to dock the ligand. A close observation revealed that Trp861 and Trp807
Fig. 9 Intermolecular interactions between iE-DAP and zNOD1-LRR in complex III. (A) Interaction of iE-DAP with the active site amino acids of zNOD1LRR shown in PyMOL. The protein, ligand and hydrogen bonds are shown as cartoon, stick and dotted lines, respectively, (B) 2-dimensional representation of iE-DAP interaction generated using the LigPlot+ program. The protein and ligand are presented in grey and green, respectively. H-bond forming amino acids are shown in grey and hydrophobic contacts as semicircles, (C) the alignment shows the key amino acid residues responsible for iE-DAP and NOD1-LLR interaction in human (red), rohu (green) and zebrafish (blue).
2950 | Mol. BioSyst., 2014, 10, 2942--2953
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
Molecular BioSystems
along with Ile863, Leu832 and Ile859 created a strong hydrophobic environment to anchor the ligand whereas polar residues such as Asp803 and Ser833 hold the ligand through high-affinity H-bond network. The other binding site residues, including Thr775, Lys777, Asn831 and Ser835 contribute crucial electrostatic interactions to the ligand. However, the 2D diagram obtained from LigPlot+67 indicates that the side chain amino group of Asn831 and the side chain protonated nitrogen atom of Trp861 interacted with iE-DAP through H-bonds (Fig. 9B). Overall, it is most likely that the loss of key H-bonds between iE-DAP and zNOD1-LRR sometime during the MD simulation is compensated by electrostatic interactions. Binding free energy calculations The energetic stability of the three trajectories of complex III was further validated by calculating the binding free energies of zNOD1-LRR and iE-DAP complexes using the MM/PBSA method. The binding of zNOD1 and iE-DAP is mostly favored by the coulombic term and is approximately equal in all trajectories of complex III. However, the second and third most contributors (nonpolar term and van der Waals energy) differed significantly. There was an approximation of 1/4th decrease in these contributions in trajectory II and III leading to comparatively stable conformations. Altogether it was found that trajectory I possessed the highest affinity between the binding partners. Therefore, trajectory I was further explored to identify the important residues responsible for the ligand binding. Computational alanine scanning revealed that most of the polar residues interacted with iE-DAP (Thr775, Lys777, Asp803 and Ser833) are crucial for high affinity ligand binding (Table 5). Among the hydrophobic residues, Ile859 was found to be essential for iE-DAP. Trp807 and Trp861, which were engaged in ligand interaction during MD simulation, presented a very negligible effect on the binding free energy. Taken together, it could be stated that zNOD1-LRR mainly recruits polar and negatively charged residues for iE-DAP recognition.
An earlier study showed that the central C-terminal region of the hNOD1 and hNOD2-LRR appeared to be a ‘hot spot’ in terms of muropeptide sensing and this work was well supported by domain swapping and site directed mutagenesis experiments.11 Interestingly, in our previous study it was proposed that LRRs 4–8 of zNOD2-LRR, located centrally to the C-terminal end is responsible for MDP recognition.40 In addition, Girardin et al. in their study showed that His778, Lys790, Gly792, Glu816, Gly818, Trp820 and Trp874 at LRRs 5–8 are indispensable for iE-DAP recognition in hNOD1.11 Similarly a study by Boyle et al. had revealed that the conserved residues on the concave surface of LRR region of hNOD1 and hNOD2 follow a common mechanism of ligand recognition.68 To map the key residues involved in zNOD1–iE-DAP interaction in our study, NOD1 sequences from zebrafish, rohu and human were aligned. As evident from Fig. 9C, the key residues which aid in ligand recognition in both human and zebrafish were conserved and differed in rohu. The reason behind the difference in iE-DAP interacting residues in rohu is not clear and warrants further investigation. In zNOD1, the conserved iE-DAP interacting residues are solvent exposed and located at the concave face of zNOD1-LRR (LRRs 4–7). Among the iE-DAP interacting residues, Lys777, Gly805, Trp807 and Trp861 were conserved in human and zebrafish, while other residues varied. Apart from the conserved residues, three additional residues Asn831, Ser833 and Ile859 were also involved in ligand recognition in zNOD1. In addition, few exceptions were also noticed, for instance, the conserved Gly779 in zebrafish corresponding to Gly792 in human was not found to interact with iE-DAP. Furthermore, Asn831, Ser833 and Ile859 of zebrafish interacted with iE-DAP, but the corresponding residues in human were excluded from the iE-DAP binding pocket. Altogether, our computational approach explained the molecular determinants responsible for ligand recognition in zNOD1. These findings suggested that the conserved LRR domains are highly specific to particular PAMPs in
Table 5 Binding free energy (kJ mol1) calculation of zNOD1-LLR and iE-DAP complexes using the MM/PBSA method: (a) wild zNOD1-LRR complex and (b) mutant zNOD1-LRR complex
Polar contribution
Nonpolar contribution
DGbind
DGcoul
DGps
DGpolar
DGvdW
DGnps
DGnonpolar
(a) Trajectory I Trajectory II Trajectory III
668.00 (5.54) 482.84 (6.03) 364.01 (5.83)
1075.75 (8.54) 1030.09 (7.42) 866.14 (7.41)
618.32 (5.94) 614.88 (5.91) 578.02 (4.72)
457.43 415.21 288.12
197.91 (3.64) 55.53 (4.73) 62.78 (4.43)
12.67 (0.04) 12.10 (0.04) 13.09 (0.03)
210.58 67.63 75.89
(b) Thr775-Ala775 Lys777-Ala777 Asp803-Ala803 Trp807-Ala807 Asn831-Ala831 Ser833-Ala833 Ile859-Ala859 Trp861-Ala861
662.54 665.99 666.82 683.52 675.37 666.45 663.73 680.34
1048.80 (8.55) 1005.69 (8.73) 954.26 (8.26) 1071.53 (8.49) 1079.76 (8.52) 1048.47 (8.49) 1074.79 (8.54) 1048.60 (8.62)
595.20 550.08 498.81 595.54 627.90 593.93 616.96 574.33
453.60 455.61 455.46 475.99 451.86 454.54 457.83 474.28
196.29 197.92 198.87 195.57 198.12 199.16 193.42 206.07
12.66 12.47 12.50 11.95 12.70 12.74 12.49 11.72
208.95 210.39 211.37 207.53 210.82 211.90 205.90 206.07
(5.45) (5.52) (5.34) (5.40) (5.47) (4.43) (5.53) (7.17)
(5.97) (5.88) (5.61) (5.91) (5.91) (5.89) (5.94) (6.10)
(3.71) (3.78) (3.63) (3.45) (3.66) (3.73) (3.71) (3.42)
(0.05) (0.05) (0.04) (0.04) (0.04) (0.05) (0.04) (0.04)
DGbind = binding free energy, DGcoul = electrostatic energy (the coulombic term), DGps = polar solvation energy, DGpolar = polar term (DGcoul + DGps), DGvdW = van der Waals energy, DGnps = nonpolar solvation energy, DGnonpolar = nonpolar term (DGvdW + DGnps), numbers in parentheses indicate standard errors. Experimental binding energy is not available.
This journal is © The Royal Society of Chemistry 2014
Mol. BioSyst., 2014, 10, 2942--2953 | 2951
View Article Online
Molecular BioSystems
different organisms and the ligand binding residues may vary from one species to another.
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Conclusion This study illustrated the modeling of zNOD1-LRR using advanced protein modeling techniques followed by long-range MD simulations. Molecular docking of iE-DAP and zNOD1-LRR followed by MD simulations of the complex indicated that the C-terminal b-sheet region (LRRs 4–7) was responsible for iE-DAP recognition, which agrees with experimental evidences in human. A combination of charged and hydrophobic residues created a compatible binding pocket on zNOD1-LRR for iE-DAP recognition. The residues His775, Lys777, Asp803, Gly805, Trp807, Asn831, Ser833, Ile859 and Trp861 were found to form the binding pocket for iE-DAP recognition. Most of the iE-DAP active site residues in NOD1 were conserved among human and zebrafish, and suggested a strong evolutionary lineage starting from lower to higher eukaryotes. The binding free energy estimation through the MM/PBSA method highlighted the importance of residues His775, Lys777, Asp803, Ser833 and Ile859 for an effective and energetically stable iE-DAP recognition. Consequently, this study unraveled mechanism of iE-DAP recognition in zebrafish, thus providing a benchmark for further investigation on NOD1-based signaling in other vertebrates with a long-term objective of finding novel therapeutic agents to conquer NOD1 related disorders.
Acknowledgements The authors thank the Director, Central Inland Fisheries Research Institute, Barrackpore, Kolkata, India, for providing institutional facility. We also thank Mr Sukanta Kumar Pradhan, Head, Department of Bioinformatics, OUAT, Bhubaneswar, and Dr Sachinandan De, Principal Scientist, NDRI, Karnal, for rendering computational facility and technical advice for improving the manuscript. JM thanks Mr Asim Kumar Jana, Senior Technical Assistant, CIFRI, Barrackpore, Kolkata, India.
References 1 R. Medzhitov and C. A. Janeway, Jr., Semin. Immunol., 1998, 10, 351–353. 2 S. Akira, S. Uematsu and O. Takeuchi, Cell, 2006, 124, 783–801. 3 E. V. Koonin and L. Aravind, Trends Biochem. Sci., 2000, 25, 223–224. ˜ ez, Nat. Rev. Immunol., 2003, 3, 4 N. Inohara and G. Nun 371–382. 5 N. Matsushima, T. Tanaka, P. Enkhbayar, T. Mikami, M. Taga, K. Yamada and Y. Kuroki, BMC Genomics, 2007, 8, 124. 6 K. Hofmann, P. Bucher and J. Tschopp, Trends Biochem. Sci., 1997, 22, 155–156.
2952 | Mol. BioSyst., 2014, 10, 2942--2953
Paper
7 A. M. Verhagen, E. J. Coulson and D. L. Vaux, Genome Biol., 2001, 2, 1–10. 8 J. Bertin and P. S. DiStefano, Cell Death Differ., 2000, 7, 1273–1274. 9 F. Martinon, K. Hofmann and J. Tschopp, Curr. Biol., 2001, 11, 118–120. 10 J. P. Ting, S. B. Willingham and D. T. Bergstralh, Nat. Rev. Immunol., 2008, 8, 372–379. ´hanno, D. Mengin-Lecreulx, P. J. Sansonetti, 11 S. E. Girardin, M. Je P. M. Alzari and D. J. Philpott, J. Biol. Chem., 2005, 280, 38648–38656. 12 M. Chamaillard, S. E. Girardin, J. Viala and D. J. Philpott, Cell. Microbiol., 2003, 5, 581–592. 13 H. Laroui, Y. Yan, Y. Narui, S. A. Ingersoll, S. Ayyadurai, M. A. Charania, F. Zhou, B. Wang, K. Salaita, S. V. Sitaraman and D. Merlin, J. Biol. Chem., 2011, 286, 31003–31013. 14 H. H. Park, Y. C. Lo, S. C. Lin, L. Wang, J. K. Yang and H. Wu, Annu. Rev. Immunol., 2007, 25, 561–586. 15 M. Tohno, W. Ueda, Y. Azuma, T. Shimazu, S. Katoh, J. M. Wang, H. Aso, H. Takada, Y. Kawai, T. Saito and H. Kitazawa, Mol. Immunol., 2008, 45, 194–203. 16 S. H. Oehlers, M. V. Flores, C. J. Hall, S. Swift, K. E. Crosier and P. S. Crosier, Dis. Models & Mech., 2011, 4, 832–841. 17 Z. Sha, J. W. Abernathy, S. Wang, P. Li, H. Kucuktas, H. Liu, E. Peatman and Z. Liu, Dev. Comp. Immunol., 2009, 33, 991–999. 18 W. Q. Chen, Q. Q. Xu, M. X. Chang, P. Nie and K. M. Peng, Fish Shellfish Immunol., 2010, 28, 18–29. 19 B. Swain, M. Basu and M. Samanta, Fish Shellfish Immunol., 2012, 32, 899–908. 20 B. R. Sahoo, B. Swain, M. R. Dikhit, M. Basu, A. Bej, P. Jayasankar and M. Samanta, Appl. Biochem. Biotechnol., 2013, 170, 1282–1309. 21 N. D. Meeker and N. S. Trede, Dev. Comp. Immunol., 2008, 32, 745–757. 22 A. Marchler-Bauer, S. Lu, J. B. Anderson, F. Chitsaz, M. K. Derbyshire, C. DeWeese-Scott, J. H. Fong, L. Y. Geer, R. C. Geer, N. R. Gonzales, M. Gwadz, D. I. Hurwitz, J. D. Jackson, Z. Ke, C. J. Lanczycki, F. Lu, G. H. Marchler, M. Mullokandov, M. V. Omelchenko, C. L. Robertson, J. S. Song, N. Thanki, R. A. Yamashita, D. Zhang, N. Zhang, C. Zheng and S. H. Bryant, Nucleic Acids Res., 2011, 39, D225–D229. 23 E. Quevillon, V. Silventoinen, S. Pillai, N. Harte, N. Mulder, R. Apweiler and R. Lopez, Nucleic Acids Res., 2005, 33, W116–W120. 24 I. Letunic, T. Doerks and P. Bork, Nucleic Acids Res., 2012, 40, D302–D305. 25 R. D. Finn, J. Mistry, J. Tate, P. Coggill, A. Heger, J. E. Pollington, O. L. Gavin, P. Gunasekaran, G. Ceric, K. Forslund, L. Holm, E. L. Sonnhammer, S. R. Eddy and A. Bateman, Nucleic Acids Res., 2010, 38, D211–D222. 26 K. Katoh and D. M. Standley, Mol. Biol. Evol., 2013, 30, 772–780. ¨ffer, R. Agarwala, S. F. Altschul, 27 G. M. Boratyn, A. A. Scha D. J. Lipman and T. L. Madden, Biol. Direct, 2012, 7, 12.
This journal is © The Royal Society of Chemistry 2014
View Article Online
Published on 29 July 2014. Downloaded by McGill University on 28/10/2014 01:28:36.
Paper
28 M. A. Kurowski and J. M. Bujnicki, Nucleic Acids Res., 2003, 31, 3305–3307. 29 O. Fornes, R. Aragues, J. Espadaler, M. A. Marti-Renom, A. Sali and B. Oliva, Bioinformatics, 2009, 25, 1506–1512. 30 Y. Yang, E. Faraggi, H. Zhao and Y. Zhou, Bioinformatics, 2011, 27, 2076–2082. 31 N. Eswar, B. Webb, M. A. Marti-Renom, M. S. Madhusudhan, D. Eramian, M. Y. Shen, U. Pieper and A. Sali, Curr. Protoc. Protein Sci., 2007, ch. 2, Unit 2.9. 32 T. Schwede, J. Kopp, N. Guex and M. C. Peitsch, Nucleic Acids Res., 2003, 31, 3381–3385. 33 C. C. Chen, J. K. Hwang and J. M. Yang, BMC Bioinf., 2009, 10, 366. 34 N. Fernandez-Fuentes, C. J. Madrid-Aliste, B. K. Rai, J. E. Fajardo and A. Fiser, Nucleic Acids Res., 2007, 35, W363–W368. ´onard, X. De Bolle and E. Depiereux, 35 C. Lambert, N. Le Bioinformatics, 2002, 18, 1250–1256. 36 J. Peng and J. Xu, Proteins, 2011, 79, 161–171. 37 S. Wu and Y. Zhang, Nucleic Acids Res., 2007, 35, 3375–3382. 38 Y. Zhang, BMC Bioinf., 2008, 9, 40. ´ll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, 39 S. Pronk, S. Pa M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854. 40 J. Maharana, M. C. Patra, B. C. De, B. R. Sahoo, B. K. Behera, S. De and S. K. Pradhan, J. Mol. Recognit., 2014, 27(5), 260–275. 41 T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092. 42 M. Wiederstein and M. J. Sippl, Nucleic Acids Res., 2007, 35, W407–W410. 43 B. Wallner and A. Elofsson, Protein Sci., 2003, 12, 1073–1086. 44 V. B. Chen, W. B. Arendall, 3rd., J. J. Headd, D. A. Keedy, R. M. Immormino, G. J. Kapral, L. W. Murray, J. S. Richardson and D. C. Richardson, Acta. Crystallogr., Sect. D: Biol. Crystallogr., 2010, 66, 12–21. 45 L. Willard, A. Ranjan, H. Zhang, H. Monzavi, R. F. Boyko, B. D. Sykes and D. S. Wishart, Nucleic Acids Res., 2003, 31, 3316–3319. 46 M. Berjanskii, P. Tang, J. Liang, J. A. Cruz, J. Zhou, Y. Zhou, E. Bassett, C. MacDonell, P. Lu, G. Lin and D. S. Wishart, Nucleic Acids Res., 2009, 37, W670–W677. 47 M. Berjanskii, Y. Liang, J. Zhou, P. Tang, P. Stothard, Y. Zhou, J. Cruz, C. MacDonell, G. Lin, P. Lu and D. S. Wishart, Nucleic Acids Res., 2010, 38, W633–W640. 48 W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
This journal is © The Royal Society of Chemistry 2014
Molecular BioSystems
49 A. K. Malde, L. Zuo, M. Breeze, M. Stroet, D. Poger, P. C. Nair, C. Oostenbrink and A. E. Mark, J. Chem. Theory Comput., 2011, 7, 4026–4037. 50 G. M. Morris, R. Huey, W. Lindstrom, M. F. Sanner, R. K. Belew, D. S. Goodsell and A. J. Olson, J. Comput. Chem., 2009, 30, 2785–2791. 51 D. Spiliotopoulos, A. Spitaleri and G. Musco, PLoS One, 2012, 7, e46902. 52 I. Massova and P. A. Kollman, J. Am. Chem. Soc., 1999, 121, 8133–8143. 53 N. A. Baker, D. Sept, S. Joseph, M. J. Holst and J. A. McCammon, Proc. Natl. Acad. Sci. U. S. A., 2001, 98, 10037–10041. 54 S. P. Brown and S. W. Muchmore, J. Med. Chem., 2009, 52, 3159–3165. 55 A. V. Kajava and B. Kobe, Protein Sci., 2002, 11, 1082–1090. 56 J. Maharana, B. Swain, B. R. Sahoo, M. R. Dikhit, M. Basu, A. S. Mahapatra, P. Jayasankar and M. Samanta, Fish Physiol. Biochem., 2012, 39, 1007–1023. 57 B. R. Sahoo, J. Maharana, G. K. Bhoi, S. K. Lenka, M. C. Patra, M. R. Dikhit, P. K. Dubey, S. K. Pradhan and B. K. Behera, Mol. Biosyst., 2014, 10, 1104–1116. 58 D. Baker and A. Sali, Science, 2001, 294, 93–96. 59 B. Kobe and J. Deisenhofer, Nature, 1995, 374, 183–186. 60 R. J. Johnson, J. G. McCoy, C. A. Bingman, G. N. Phillips Jr. and R. T. Raines, J. Mol. Biol., 2007, 368, 434–449. 61 D. W. Buchan, F. Minneci, T. C. Nugent, K. Bryson and D. T. Jones, Nucleic Acids Res., 2013, 41, W349–W357. 62 G. N. Ramachandran, C. Ramakrishnan and V. Sasisekharan, J. Mol. Biol., 1963, 7, 95–99. 63 R. A. Laskowski, M. W. MacArthur, D. S. Moss and J. M. Thornton, J. Appl. Crystallogr., 1993, 26, 283–291. ¨thy, J. U. Bowie and D. Eisenberg, Nature, 1992, 356, 64 R. Lu 83–85. 65 C. Colovos and T. O. Yeates, Protein Sci., 1993, 2, 1511–1519. 66 T. Tanabe, M. Chamaillard, Y. Ogura, L. Zhu, S. Qiu, J. Masumoto, P. Ghosh, A. Moran, M. M. Predergast, ˜ ez, ´n G. Tromp, C. J. Williams, N. Inohara and G. Nu EMBO J., 2004, 23, 1587–1597. 67 R. A. Laskowski and M. B. Swindells, LigPlot+: multiple ligand-protein interaction diagrams for drug discovery, J. Chem. Inf. Model., 2011, 51, 2778–2786. 68 J. P. Boyle, S. Mayle, R. Parkhouse and T. P. Monie, Front. Microb. Immunol., 2013, 4, 317.
Mol. BioSyst., 2014, 10, 2942--2953 | 2953
